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(ABSTRACT) 

Because of the inherent complexity of fiber-reinforced laminated composites, it can 
be challenging to manufacture composite structures according to their exact design 
specifications, resulting in unwanted material and geometric uncertainties. In this 
research, we focus on the deterministic and probabilistic stability analysis of lam- 
inated structures subject to subtangential loading, a combination of conservative 
and nonconservative tangential loads, using the dynamic criterion. 

Thus a shear-deformable laminated beam element, including warping effects, is 
derived to study the deterministic and probabilistic response of laminated beams. 
This twenty-one degrees of freedom element can be used for solving both static and 
dynamic problems. In the first-order shear deformable model used here we have 
employed a more accurate method to obtain the transverse shear correction factor. 
The dynamic version of the principle of virtual work for laminated composites is 
expressed in its nondimensional form and the element tangent stiffness and mass 
matrices are obtained using analytical integration. The stability is studied by giving 
the structure a small disturbance about an equilibrium configuration, and observing 
if the resulting response remains small. 

In order to study the dynamic behavior by including uncertainties into the 
problem, three models were developed: Exact Monte Carlo Simulation, Sensitivity- 
Based Monte Carlo Simulation, and Probabilistic FEA. These methods were inte- 
grated into the developed finite element analysis. Also, perturbation and sensitivity 
analysis have been used to study nonconservative problems, as well as to study the 
stability analysis using the dynamic criterion. 
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Chapter 1 

Analysis of Laminated Structures 
with Uncertain Properties 


The increased use of fiber-reinforced laminated composites has created a new 
interest in the analysis of laminated structural elements such as laminated beams, 
plates, and shells. A better understanding of the structural behavior of these com- 
posite structures is much needed. Mechanical and physical properties of laminated 
composite structures can become uncertain due to changes in various factors like 
fiber orientations, curing temperature, pressure and time, voids, and impurities 
among others. Thus, for this kind of material it is extremely important to identify 
and use the most appropriate model of uncertainty. 

Although researchers in the past have studied the kinetics of fiber-reinforced 
laminated composites, in the present work a new twenty-one degree of freedom beam 
element is derived to perform static, dynamic, and stability analysis through the 
finite element nonlinear analysis. This encompasses the first goal of this research. 

The second goal is the probabilistic study of the structural behavior by includ- 
ing uncertainties into the problem through the probabilistic finite element method 
(PFEM), exact Monte Carlo simulation (EMCS), and sensitivity- based Monte Carlo 
simulation (SBMCS). 

The purpose of this chapter is to: (i) introduce uncertainties in laminated com- 
posites, (ii) present previous work done on finite element nonlinear analysis of 
laminated beams and PFEM, and (iii) state the motivation and explain the scope 
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Figure 1.1: Fiber-reinforced laminated composites. 


of the present work. 


1.1 Analysis of Laminated Composites 

The composites considered here are fiber- reinforced laminated composites, a hy- 
brid class of composite materials involving both fibrous composite materials and 
lamination techniques. Fiber-reinforced laminated composites are of great interest 
to aerospace engineers because these composites can be tailored to best match the 
design requirements of a specific structural application, while allowing structural 
components to remain lightweight. These laminated structures are obtained b\ 
stacking two or more laminae with their fibers oriented in different directions, as 
shown in Fig. 1.1. 

The analysis of fiber-reinforced laminated composites is far more complex when 
compared to conventional materials because composites are inhomogeneous through 
the thickness and generally anisotropic. Anisotropic materials exhibit directional 
characteristics and thus bring shear-extension coupling into the analysis. When 
laminates are unsymmetrically stacked, bending-stretching coupling must be in- 
cluded in the analysis. Moreover, differences in elastic properties between fiber 
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filaments and matrix materials lead to inplane-shear coupling. These couplings 
only add further complexity to the solution. 

The high ratio of inplane modulus to transverse shear modulus makes the clas- 
sical lamination theory, which neglects the effect of out-of- plane strains, inadequate 
for the analysis of multilayer composites (Mallikarjuna and Kant, 1993). In such 
a case, one should use a theory which includes transverse shear deformation; an 
example is the first order shear deformation theory (FSDT). In fact, the presence 
of shear forces introduces shear stresses in a beam which may cause appreciable 
angular deformation of the beam. As a consequence, transverse normals no longer 
remain normal after deformation. In laminated composites, such a behavior can be 
very important. Because laminated composite materials have very low transverse 
shear modulus compared to their in-plane moduli, the classical laminated theory 
will under-predict deflections and over-predict frequencies and buckling loads (Singh 
et al., 1991). 

The FSDT, which ignores the effects of cross-sectional warping, leads to an 
unrealistic (constant) variation of the transverse shear stress through the laminate 
thickness. Higher-order shear deformable theories (HSDT) do not require the use 
of shear constants because they model in a realistic manner the parabolic variation 
of the transverse shear stress through the laminate thickness. However, for the 
type of problems considered in this research, the FSDT is used. When using the 
FSDT, the shear correction factors should be used (Reissner, 1945; Whitney, 1973; 
Reissner, 1972). Here a procedure similar to that of Cohen (1978) is used. 

St. Venant’s theory, as it concerns torsion, assumes that the cross-sections of 
the beam maintain their original shape, although they are free to warp in the 
axial direction. The warping displacement is postulated to be proportional to the 
rate of twist, which is assumed to be constant along the beam axis (Shames and 
Dym, 1985). However, these classical theories of bending and torsion may result in 
erroneous predictions, especially when the warping constraint is present. 
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Within the warping constraint beam model, the rate of twist can not be assumed 
to remain constant along the axis of the beam. Thus for an arbitrary cross section 
we must abandon the assumption that plane sections remain plane and introduce 
a warping function to account for warping in the x-direction. Many researchers 
have developed more realistic models (Librescu and Song, 1991; Fang and Springer, 
1991; Berdichevsky et al., 1992). They have shown that for anisotropic cantilevered 
beams, the warping inhibition induced by the restraint of torsion suggests that 
the St. Venant’s twist concept is not applicable. Crawley and Dugundji (1980) 
studied the warping restraint effect on torsional vibration with bending-torsion 
coupling. Later, Kaza and Kielb (1984) investigated the torsional vibration of 
rotating pretwist beams by varying the aspect ratio and including the warping 
effect. Librescu and Song (1991) incorporated the effects of primary and secondary 
warping restraint in the study of composite thin-walled beam structures. Although 
many researchers have suggested various methods of approximating the warping 
distribution, here we propose to use a simpler method to calculate the warping 
function for thin beams. 

Although the analysis of a laminated composite structure as a three-dimensional 
solid would be most accurate, it is numerically costly. Without loss of accuracy, 
thin-walled laminated composites can be modeled as two-dimensional structures 
(i.e., plates and shells) because they (i) are made with their planar dimensions one 
or two orders of magnitude larger than their thickness, and (ii) are often used in 
applications that require only axial and bending strength. Moreover, structures 
having high length-to-width and length-to-thickness ratios can be approximated 
using one-dimensional theories, as is the case of helicopter rotor blades and some 


wings. 
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1.2 Finite Element Nonlinear Analysis 


The study of fiber-reinforced laminated composite structures using discretized meth- 
ods such as the finite element analysis is of great interest to researchers. With few 
exceptions, most researchers have basically considered the traditional approach of 
the finite element analysis, one of using numerical integration. An advantage of 
numerical integration is that it is easier to implement in an existent finite element 
code. On the other hand, its use is disadvantageous because it may lead to shear- 
locking and in some cases to convergence problems (Bathe, 1996). 

Only a few researchers have tried to integrate symbolic computation into the 
finite element models (Yang, 1994; Teboub and Hajela, 1995). An advantage of 
symbolic integration is that it saves computational time because the tangent stiff- 
ness matrix and the internal force vector, per iteration, for each element are eval- 
uated only once. To the best of our knowledge, no work that integrates symbolic 
manipulators into the finite element nonlinear analysis for laminated composites is 
available. 

In this research, it is of interest to study structures that undergo large defor- 
mations. Geometric nonlinearities must be included when large deformations are 
important. Moreover, for theories including shear deformation and geometric non- 
linearities, a two-dimensional theory would require a large number of degrees of 
freedom. A one-dimensional model is desirable and thus used here. 

In the last two decades, some attention has been paid to the development of lam- 
inated beam elements. Earlier, a 12-dof element was developed and formulated for 
deterministic symmetric laminated beams to study their static and dynamic behav- 
iors (Chen and Yang, 1985). Later a 20-dof element (Kapania-Raciti Element) was 
developed to study static, free vibration, buckling, and nonlinear vibrational anal- 
ysis of unsymmetrically laminated beams (Kapania and Raciti, 1989a). Kapama 
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and Raciti (1989b, 1989c) gave an extensive review of the literature on laminated 
beams and plates that existed prior to 1989. However, no derived beam elements 
available before or since 1989 have all of the following features in one element: 
use of symbolic integration, all three displacements (axial, transverse, and lateral), 
torsional and warping effects, inplane rotation, and shear effects. 

Bassiouni et al. (1999) used a finite element model to obtain the natural frequen- 
cies and mode shapes of laminated composite beams. The finite element consisted 
of five nodes: the two end points, one at one-third, one at two-thirds, and the mid- 
point. The displacement components are lateral displacement, axial displacement, 
and rotational displacement. Closed form solutions for stiffness and mass matri- 
ces are given. Although this is one of the most complete beam elements, it does 
not take into account inplane shear, shear deformations, transverse deflections, and 
torsional effects. 

Finite elements based on the FSDT are challenging because they add additional 
degrees of freedom to the problem. However, quite a few researchers have developed 
such elements. Koo and Kwak (1994) proposed a finite element suitable for the 
analysis of composite frames based on the first-order shear deformation theory. 
The deflection is separately interpolated for bending and shearing with cubic and 
linear functions, respectively. Carrera and Villani (1994) dealt with the nonlinear 
analysis of multilayered axially compressed plates in static conservative cases. 

A formulation for the exact dynamic stiffness matrix for symmetric and unsym- 
metrically laminated beams has been derived using the exact shape functions for the 
deflection and bending slope of composite laminated beam elements (Abramovich 
et al., 1995; Eisenberger et al., 1995). The formulation is based on the FSDT and 
includes rotatory inertia effects developed by Abramovich (1992). 
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Some researchers have analyzed laminated beams using the state-space solu- 
tion. Teboub and Hajela (1995) studied the free vibrations of anisotropic lami- 
nated composite beams using FSDT, including inplane and rotatory inertias, using 
a state-space solution instead of the finite element method. Later, Khdeir and 
Reddy (1997) used the state-space concept in conjunction with the Jordan canon- 
ical form to solve the governing equations for the bending of cross-ply laminated 
beams. They used classical, second-order, and third order theories to develop exact 
solutions for symmetric and unsymmetric cross-ply laminated beams. The solution 
is mainly based on the theories previously employed to investigate the vibration and 
buckling analysis of cross-ply laminated beam (Khdeir and Reddy, 1994; Khdeir, 
1996). 

Some research has been focused on deriving beam elements for laminated com- 
posites using higher order shear deformation theories (HSDT). Manjunatha and 
Kant (1993) developed a set of higher-order theories for the analysis of composite 
and sandwich beams by using C° finite elements. By incorporating a more realistic 
nonlinear displacement variation through the beam thickness, they eliminated the 
need for shear correction coefficients. Kam and Chang (1992) studied the bending 
and free vibration behavior of laminated composite beams using FSDT and HSDT. 

Shi et al. (1998) investigated the influence of the interpolation order of the ele- 
ment bending strain on the solution accuracy of composite beam elements derived 
using HSDT and presented a simple and accurate third-order composite beam el- 
ement. They concluded that the strain expression that gives the higher order of 
bending strain interpolation should be chosen for the finite element modeling of 
composite beams based on HSDT. 

Kadivar and Mohebpour (1998) studied the finite element dynamic response of 
an unsymmetrically laminated composite beam subject to moving loads. The one- 
dimensional finite element is derived based on classical lamination theory, first-order 
shear deformation theory, and higher-order shear deformation theory. 
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Subramanian (2001) developed a two-noded C 1 finite element of eight degrees 
of freedom per node, using a HSDT, for flexural analysis of symmetric laminated 
composite beams. Lam and Zou (2001) developed a higher-order shear deformable 
finite strip for the analysis of composite laminates. The formulation allows C 
continuity with nine variables and can be used to analyze both symmetrically and 
unsymmetrically laminated plates. 

Yildirim and Kiran (2000) studied the out-of-plane free vibration problem of 
symmetric cross-ply laminated beams by the transfer matrix method. The formu- 
lation is based on the first-order shear deformation theory. In their formulation it is 
possible to isolate the effects of rotatory inertia, transverse shear, and axial defor- 
mations to study their influence on the natural frequencies. In a subsequent work, 
Yildirim (2000) used the stiffness method for the solution of the purely inplane 
free vibration problem of symmetric cross-ply laminated beams. In the former, 
six degrees of freedom were defined for an element, four displacements and two 
rotations. 

Others have used layerwise theories in the development of laminated beam ele- 
ments. Averill and Yip (1996) developed an accurate, simple, and robust two-noded 
C° finite elements based on shear deformable and layerwise (zig-zag) laminated 
beam theories. The two-noded element has only four degrees of freedom per node. 
The formulation is only valid for small deformations. 

Cho and Averill (1997) developed a beam finite element based on a new dis- 
crete layer laminated beam theory with sublaminate first-order zig-zag kinematic 
assumptions for both thin and thick laminated beams. The finite element is de- 
veloped with the topology of a four-noded rectangle, allowing the thickness of the 
beam to be discretized into several elements, or sublaminates. 

Only few researchers have worked on nonlinear analysis of laminated beams 
using the finite element method. Murin (1995) formulated a nonlinear stiffness 


1.2. FINITE ELEMENT NONLINEAR ANALYSIS 


9 


matrix of a finite element without making any simplifications. The matrix includes 
the quadratic and cubic dependencies of the unknown increments of the generalized 
nodal displacements into the initially linearized system of equations. However, the 
formulation is limited to isotropic materials and not applied to laminated compos- 
ites. 

Patel et al. (1999) studied the nonlinear flexural vibrations and post-buckling of 
laminated orthotropic beams resting on a class of two parameter elastic foundations 
using a three-noded shear flexible beam element. 

After reviewing most of the work done in this field, a number of challenges 
intrigued us. In some formulations, the transverse displacement was split as the 
addition of shear and bending contributions. This leads to a singular mass matrix 
when considering rotatory inertias. On the other hand, one can also consider the 
total transverse deflection and introduce shear using the rotation of the normal to 
the middle surface in the displacement field. 

Another fact is that when considering in- plane shear rotation 0 = dU/dy, the 
lateral displacement v cannot be ignored. Ignoring the lateral displacement brings 
inconsistencies in the strain-displacement relations. 

Moreover, when studying the uncertainties of symmetrically and unsymmetri- 
cally laminated beams, one should have a beam element capable of varying material 
and geometric properties. The present element enables us to vary all the properties 
of the beam, a characteristic which has helped us to study the probabilistic nature 
of laminated beams in the presence of uncertainties (Kapania and Goyal, 2001, 
2002 ). 

To the best of the author’s knowledge, there is no shear deformable laminated 
beam element for the study of large deflection, including torsional and warping 
effects, lateral displacement, and inplane displacement as a degree of freedom. 
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Thus, there seems to be a need for a much improved laminated beam element 
for the static and dynamic analysis of large deflections of symmetric and unsym- 
metrically laminated composite beams. An element that does not require the use of 
numerical integration but uses symbolic manipulators such as Mathematica would 
be most helpful in a study of composite beams with uncertainties. 


1.3 Dynamic Stability of Laminated Beams 

Usually structures may become unstable under an increasing compressive load. For 
such problems, stability analysis becomes a significant tool when it comes to design 
and analysis of structures such as aircraft, ships, and automobiles. 

The stability of the equilibrium state can be studied using three criteria: (i) 
static or adjacent equilibrium criterion, (ii) energy criterion, (iii) dynamic (or ki- 
netic or vibration) criterion. For a conservative system, the only possible initial 
instability is of divergence type. For a nonconservative system, however, instability 
can be of divergence, flutter, or both, depending on the amount of nonconserva- 
tiveness. Since only the kinetic method works for both conservative and nonconser- 
vative systems, it will be the method used in this research. The dynamic criterion 
relates the critical loads and the natural frequencies of the system. 

Stability analysis using the dynamic criterion has been used by various re- 
searchers. Argyris and Symeondis (1981) presented a nonlinear finite element anal- 
ysis of elastic structures subject to nonconservative forces. They derived a general 
theory to study the stability behavior of non-self-adjoint boundary-value problems. 

Hasegawa et al. (1988) studied the elastic instability and the nonlinear finite dis- 
placement behavior of spatial thin-walled members under displacement-dependent 
loads. They presented a general formulation to derive the loading stiffness matrix. 
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The dynamic stability of structures under a follower force using the finite ele- 
ment method has also been studied by Chen and Yang (1989), Chen and Ku (1991a, 
1991b), Saje and Jelenic (1994), Vitaliani et al (1997), Kim and Kim (2000), and 
Detinko (2001). 

Moreover, there has been some interest to study structures subject to both 
conservative and nonconservative loads. When the nonconservativeness is added as 
a fraction of the purely tangential follower loads, it is called a subtangential load. 
Rao and Rao (1987a, 1987b) studied the stability of a cantilevered beam under 
a subtangential follower load using the static and dynamic criteria. Gasparini 
et al. (1995) discussed the transition between the stability and instability of a 
cantilever beam subjected to a partially follower load by using the FEM. 

Later, Zuo and Schreyer (1996) studied the instability of a cantilevered beam 
and a simply supported plate, subjected to a combination of fixed and follower 
forces. They introduced a nonconservative parameter to account for all possible 
combinations of these forces. They showed that for the beam, instability changes 
from divergence to flutter at a critical value of this parameter; for values of the 
parameter above the critical value, the flutter instability remains as the only insta- 
bility pattern; and for the plate, the instability is governed by flutter for a certain 
range of the nonconservative parameter, even though divergence instability still 

exists. 

Ryu et al. (1998) investigated the dynamic stability of cantilevered Timoshenko 
vertical columns having a tip rigid body and subjected to the action of subtangential 
forces. They refer to subtangential force as a combination of tangential follower 
force with the vertical force. 

Only few researchers have studied the stability of laminated structures under 
a tangential follower loads, such as Xiong and Wang (1987). They presented an 
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analytical method for calculating the stability of a laminated column under a Beck- 
type load including shear deformation and rotatory inertia. 

In fact, aerospace structures may not only be subject to conservative loads but 
also to nonconservative loads as well. To the best of the author s knowledge, no 
work has been found on the stability of laminated beams subject to subtangential 
loading using the dynamic criterion. Thus, we will study how the stability of 
laminated composite beams is affected by subtangential loads. 


1.4 A Computational Probabilistic Analysis 


Because of the inherent complexity of composite materials, fiber-reinforced lam- 
inated structures can be difficult to manufacture according to their exact design 
specifications, resulting in unwanted uncertainties. In fact, during the manufactur- 
ing of laminates, material defects such as interlaminar voids, delamination, incorrect 
orientation, damaged fibers, and variation in thickness may be introduced (Reddy, 
1997). 

The design and analysis using conventional materials is easier than those using 
composites because for conventional materials both material and geometric prop- 
erties have either little or well known variation from their nominal value. On the 
other hand, the same cannot be said for the design of structures using laminated 
composite materials. Thus, the understanding of uncertainties in laminated struc 
tures is highly important for an accurate design and analysis of aerospace and other 
structures. Elishakoff (1998) has suggested three different approaches to study un- 
certainties: (i) probabilistic methods, (ii) fuzzy set or possibility-based methods, 
and (iii) anti-optimization. 

The non-probabilistic methods such as fuzzy set theory and anti-optimization 
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are used when data regarding the uncertain parameter are not available or little 
is known about the probability density function (Ayyub, 1994; Elishakoff, 1995). 
However, these uncertainties are ignored in this research because non-probabilistic 
methods are beyond the scope of the present work. 

In this research, the noncognitive sources of uncertainty are of great interest 
and are treated using probabilistic methods or non-probabilistic methods. The 
noncognitive sources of uncertainty (i.e., material and geometric variations) are in 
general quantified and information about the uncertainty of these parameters may 
be available. When sufficient data are available to predict the probability density 
function, then a probabilistic method can be used. Thus, throughout this disser- 
tation uncertainties due to noncognitive sources are studied using a probabilistic 
approach. 

The randomness of noncognitive sources leads to variations in the stiffness and 
mass coefficients of the laminates. These uncertainties may involve geometric quan- 
tities (e.g., ply-orientations and dimensions), material properties (e.g., elastic mod- 
ulus, shear modulus, Poisson’s ratio, and material density), and external properties 
(e.g., thermal and loading effects). However, in this research only those uncertain- 
ties involving material and geometric properties are considered. 

The probabilistic analysis can be performed using either an analytical or a com- 
putational approach. An analytical approach would be most accurate although 
cumbersome and impractical except for very simple systems. However, with the 
availability of high-speed computers, the finite element method has become a stan- 
dard tool for engineers to analyze structures with complex geometry, including 
various sources of nonlinearities. However, the deterministic finite element method 
fails to take into account uncertainty in different parameters of the structure, and 
thus cannot be used for reliability analysis. 
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Various methods exist to analyze an uncertain structure by integrating prob- 
abilistic aspects into the finite element modeling. Especially, there has been a 
growing interest in applying these methods to better understand laminated com- 
posite structures by integrating the stochastic nature of the structure in the finite 
element analysis (Schueller, 1997). When the probabilistic nature of material prop- 
erties, geometry, and/or loads is integrated in the finite element method, such a 
concept is called probabilistic finite element method (PFEM). 

The probabilistic finite element analysis (PFEA) can be classified into two cat- 
egories: perturbation techniques and simulation methods. Perturbation techniques 
are based on series expansion (e.g., Taylor Series) to formulate a linear or quadratic 
relationship between the randomness of the material, geometry, or load and the ran- 
domness of the response (Nakagiri and Hisada, 1988a; 1988b). Simulation methods 
such as Monte Carlo simulation rely on computers to generate random numbers 
from the material, geometry, or load uncertainties and correlate the probabilistic 
response to it (Shinozuka, 1972; Fang and Springer, 1993; Vinckenroy et al., 1995). 

A considerable amount of research has been made in the field of random struc- 
tures using the stochastic finite element method. Conteras (1980) and Vanmarcke 
et al. (1986) applied the method to the analysis of static and dynamic problems. 
Collins and Thompson (1969) applied it to the analysis of eigenvalue problems. 
Kiureghian and Ke (1988) and Zhang et al. (1996) applied perturbation methods 
to structural design. The major application has been for design purposes and all 
the work has been applied to isotropic materials. 

Chakraborty and Dey (1995) developed a stochastic FEM for the analysis of 
structures having statistical uncertainties in both material properties and exter- 
nally applied loads. These uncertainties were modeled as homogeneous Gaussian 
stochastic processes. They used the Neumann expansion technique to invert the 
stochastic stiffness matrix. They did not consider a stochastic mass matrix, and 
the formulation was applied only to linear stiffness matrices. Also, no application 
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was made to composite materials. 

The probabilistic analysis requires the derivatives of the structural matrices as 
well as the derivatives of the eigenvalues, eigenvectors, and displacements. Lee and 
Lim (1997) presented an approach for extending sensitivity methods to include the 
structural uncertainty with random parameters using perturbation techniques. 

Derivatives of eigenvectors with respect to design variables are very useful in 
certain analyses and design applications. Many researchers have developed various 
sensitivity-based methods to calculate these derivatives (Fox and Kapoor, 1968, 
Plaut and Huseyin, 1973; Haftka and Adelman, 1986; Liu et ah, 1995). The sen- 
sitivity analysis and calculation of laminated composites as a tool for design opti- 
mization has been studied by various researchers such as Pederson (1987), Mateus 
et al. (1991), and Chen et al. (1996). 

Brenner and Bucher (1995) presented a stochastic finite element- based reliabil- 
ity analysis of large nonlinear structures under dynamic loading, involving both 
structural and loading randomness, with relatively little computational effort when 
compared to the traditional Monte Carlo methods. Papadopoulos and Papadrakakis 
(1998) used a weighted integral method in conjunction with Monte Carlo simulation 
for the stochastic finite element-based reliability analysis of space frames. 

Chakraborty and Dey (1996) implemented the stochastic finite element simula- 
tion of random structures on uncertain foundations under random loading. Later, 
Chakraborty and Dey (1998) proposed a stochastic FEM for the frequency domain 
for the analysis of structural dynamic problems involving uncertain parameters. 

Recently, Oh and Librescu (1997) studied the free vibration and reliability of 
cantilever composite beams featuring structural uncertainties. They used a stochas- 
tic Rayleigh-Ritz formulation. Graham and Deodatis (2000) studied the variability 
of the response displacements and eigenvalues of structures with multiple uncertain 
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material and geometric properties. 

Imai and Frangopol (2000) reviewed the theory of finite element reliability anal- 
ysis of geometrically nonlinear elastic structures based on the total Lagrangian for- 
mulation. They also provided developments in computer implementation and es- 
tablished the basis of understanding of the applications presented in a subsequent 
work (Frangopol and Imai, 2000). 

Mei et al. (1998) used a wavelet-based stochastic analysis to analyze isotropic 
beam structures. Sobczyk et al. (1996) analyzed the dynamics of structural systems 
with randomly varying parameters using the random integral equation theory. 

The presence of noncognitive uncertainties will lead to randomness in the ma- 
terial and geometric parameters. Because of the uncertain nature of material and 
geometric parameters, the stability and vibration analysis of the trivial equilibrium 
state (assuming it exists) will be affected as well. 

The analysis of structures under random loading has been studied for a large 
class of problems (Maymon, 1998). Moreover, noncognitive uncertainties in com- 
posite material properties have been studied by Nakagiri and Hisada (1983), Nak- 
agiri et al. (1987), Ibrahim (1987), Leissa and Martin (1990), and Oh and Librescu 
(1997). 

In the dynamic analysis of the present problem, the random nature of the stiff- 
ness matrix, mass matrix, eigenvalues, and eigenvectors can be studied using a Tay- 
lor series expansion up to second order about the mean of each random variable. 
This approach was recently used by Zhang and Ellingwood (1995) to solve buck- 
ling problems. Oh and Librescu (1997) used a similar formulation for a stochastic 
Rayleigh- Ritz approach to study the free vibrations of laminated composites. 

However, to the best of the author’s knowledge, no work has been found in 
integrating the probabilistic finite element method into the vibrations and stability 
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analysis using perturbation methods. Thus, this work will intend to integrate both 
the existent methods of analyzing the stability of the equilibrium state using the dy- 
namic criterion and three different methods for probabilistic analysis: probabilistic 
finite element method, sensitivity-based Monte Carlo simulation, and Monte Carlo 
simulation. 


1.5 Present Work 

The goals of this research can be summarized as follows: 

1. Develop a shear deformable laminated beam element using the first-order 
shear deformation beam theory featuring geometric and material uncertain- 
ties. 

2. Perform the deterministic stability analysis of laminated beams subjected to 
subtangential loading using the dynamic criterion. 

3. For a structure with imperfections in the ply angles and ply axial modulus of 
elasticity, perform stability analysis of uncertain systems using the dynamic 
criterion of conservative and nonconservative systems. 

As a first step in this journey to study the uncertainties of laminated structures, 
it is assumed that the random variables are dependent on only one dimension, i.e., 
the x-axis. This assumption requires the use of a one-dimensional finite element 
model. The elements previously developed in the literature have raised our curiosity 
because of the following reasons: 


• In some cases, if the transverse displacement is split as the addition of shear 
and bending contributions, it leads to a singular mass matrix when inertia 
terms are kept. 
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• The probabilistic nature of laminated beams cannot be easily studied. 

• Most derivations use numerical integration as opposed to symbolic manipu- 
lators. 

• When considering inplane shear rotation {3 = dU/dy and twist angle r = 
dW/dy, the lateral displacement V cannot be ignored. 

Thus a shear deformable laminated beam element, with twenty-one degrees-of- 
freedom, is developed using first-order shear deformation beam theory to account for 
uncertainties. In Chapter 2, we present the Total Lagrangian formulation for gener- 
ally anisotropic laminated composite structures and state the generalized principle 
of virtual work (GPVW). We also nondimensionalize the GPVW for our problem. 

Chapter 3 is devoted to the formulation of the shear deformable beam element 
which is obtained by discretizing the generalized principle of virtual work. Since it 
is desired to study the stability of the equilibrium state using the dynamic criterion, 
the equations of motion are perturbed about a trivial equilibrium state, we present 
the finite element formulation for the stability analysis of partially conservative 
systems. In Chapter 4 we present the deterministic analysis of laminated structures 
subjected to subtangential loading. 

The second part of this research is dedicated to develop and study the stability 
of symmetrically and unsymmetrically laminated beams featuring mechanical and 
geometric uncertainties. Thus the probabilistic method using the finite element 
method is developed in Chapter 5. Results and discussions for the probabilistic 
analysis are presented in Chapter 6. 

The final chapter is a brief summary of this dissertation and includes our final 
remarks. The last section of this chapter includes areas in which this work can be 
expanded. 


Chapter 2 

Equations of Motion for Generally 

Anisotropic Laminated 

Composites 


The equations of motion can be derived using energy methods and the principle 
of virtual work. The advantage of the principle of virtual work over the principle 
of minimal total potential energy is that the formulation is applicable to both 
conservative and nonconservative problems. 

In the present work, we are interested in studying the stability of structures 
subject to conservative and nonconservative loads. Methods based on the principle 
of minimum potential energy are not helpful because they are not applicable to 
structures subjected to nonconservative loads, such as the follower loads. Moreover, 
the stability analysis of this kind of problem can only be studied using the dynamic 
criterion. Thus, here the generalized principle of virtual work which includes 
inertial loads — is preferred. 

In this chapter, the equations of motion are presented using the generalized 
principle of virtual work for the static and dynamic analysis of both symmetri- 
cally and unsymmetrically laminated beams with uncertain parameters subject to 
conservative and nonconservative loading. 
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2.1 Basic Assumptions 

Although using general theories and equations would be most accurate, they would 
be impractical for a wide range of problems. In this work, several assumptions 
simplify these theories, helping us grasp a better understanding of the problem. 

As a first simplification, thermal and piezoelectric effects are ignored. Every 
structure undergoes temperature changes. However, we will restrict ourselves to 
cases for which the temperature gradient is very small or negligible. The structures 
considered here do not have any piezoelectric devices. 

We will consider a general out-of-plane warping. Although many researchers 
have ignored warping, as the present results show, it cannot be ignored for certain 
types of laminated beams. 

Although a two-dimensional theory would be most accurate, the discretization 
would require a large number of degrees of freedom. Structures having one dimen- 
sion far larger than the other two can be approximated using one-dimensional the- 
ories. In the present work, the plate is assumed to have a high length-to-width and 
length-to-thickness ratios. Thus the plate can be modeled using a one-dimensional 
theory, such as beam theory. Here we use the term IdTuiTiuted becwi to refer to plate 
strips. 

For most aerospace applications, the thickness of the plate can be assumed far 
smaller than the length and the width. As a result, the magnitudes of the stresses 
acting on the surface parallel to the mid-plane are small compared to the bending 
and membrane stresses. As a result, the state of stress can be approximated as a 
state of plane stress. 

Laminated composite materials are made of fiber-reinforced lamina of different 
properties. It is assumed that each lamina is a continuum (i.e., no empty spaces, 
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z, w 



Figure 2.1: Reference coordinate system 

voids, internal delaminations, or material defects exist) and it behaves as a linear 
hyperelastic material. 

Fiber-reinforced lamina are often characterized by materials having three mu- 
tually orthogonal planes of material symmetry referred to the principal material 
directions (i.e., orthotropic materials). Therefore, each lamina is assumed to be 
characterized by orthotropic materials. 

In laminated plates with large bending-stretching coupling, the effect of curva- 
ture may arise during the curing process of an unsymmetrically laminated plate. 
However, this effect is ignored and the present formulation assumes a perfectly 
straight laminated plate in its initial configuration. In the case of post-buckling 
analyses, imperfections can be included by assuming the loading to be eccentric. 

A cartesian coordinate system x-y-z, shown in Fig. 2.1, is used and it is placed 
at the mid-surface of the laminate: the axial displacement u is associated with 
the x-axis which lies along the length of the beam, the lateral displacement v is 
associated with the y - axis which is placed at the middle of the beam, and the 
transverse displacement w (upwards) associated with the 2 -axis which is placed at 
the midsurface. 
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Further, it is assumed that the x-z plane divides the beam in two identical parts: 
in other words, material, geometry (with the exception of ply angles), and loading 
are symmetric about the x-z plane. 


2.2 Nonlinear Behavior 

Many structural problems undergo nonlinear behavior, even though sometimes such 
behavior may be difficult to identify. In general, there are four ways in which 
structural nonlinear behavior can occur (Doyle, 2001). These are: 

1. Material nonlinearity. The material stress-strain relationship is actively 
nonlinear. In this case, material behavior depends on the current deformation 
state and possibly past history of the deformation. Material nonlinearity can 
be observed in structures undergoing nonlinear elasticity, plasticity, nonlin- 
earity viscoelasticity, creep or other inelastic effects. 

2. Geometric nonlinearity. There is a nonlinear strain-displacement relation- 
ship. The change in geometry, as the structure deforms, is taken into account 
when forming the strain-displacement relationship and hence the equilibrium 
equations. Geometric nonlinearity may be due to large strains or small strains 
but with large displacements and/or rotations (cables, leaf-springs, arches, 
fishing rods, snap-through buckling). 

3. Application of nonlinear forces. The magnitude or direction of the 
applied forces changes with application to the structure ( nonlinear force- 
deflection relationship). This could be due to pressure loadings, gyroscopic 
forces, or follower forces. 
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4. Displacement boundary condition nonlinearity. The displacement 
boundary conditions depend on the deformation of the structure (nonlinear 
displacement-deformation relationship ). The most important application is 
in contact problems, where the displacement is highly dependent on the rela- 
tionship between two contact surfaces (normal force and friction present). 

Throughout this dissertation, nonlinearities in both material and boundary are 
neglected and only geometric and force nonlinearities are considered. We use ge- 
ometric nonlinearities to analyze laminated structures with large displacements, 
small strains, and moderate rotations. The types of force nonlinearities considered 
here are those due to tangential follower loads. 


2.3 Kinematic Description 

When geometric nonlinearities are included, three kinematic descriptions are avail- 
able: 

1. Total Lagrangian Description (TL): the reference configuration is sel- 
dom or never changed and it is often kept equal to the base configuration 
throughout the analysis. 

2. Updated Lagrangian Description (UL): strains and stresses are redefined 
as soon as the reference configuration is updated. 

3. Corotational Description (CR): strains and stresses are measured from 
the corotated configuration, whereas the base configuration is maintained as 
reference for measuring rigid body motions. 

In problems related to uncertainties, usually the noncognitive uncertainties are 
known in the structure’s reference configuration. In the deformed configuration 
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the knowledge of these uncertainties is most likely unknown. A kinematic descrip- 
tion that seldom changes its reference configuration is desirable. Thus, the Total 
Lagrangian description, a widely used formulation in continuum-based analysis, is 
used (Bathe, 1996). 


2.4 Continuum Mechanics 


Continuum mechanics is essential for nonlinear analysis. This section starts with 
the displacement field used in this dissertation. Next, we provide a brief description 
of deformation, present the Green-Lagrange strain measures, and finally express the 
Second Piola-Kirchhoff stresses (PK2) in terms of the physical (Cauchy) stresses. 

2.4.1 Displacement Field 


The analysis of fiber-reinforced laminated composites is far more complex when 
compared to conventional materials because composites are inhomogeneous through 
the thickness and generally anisotropic. Anisotropic materials exhibit directional 
characteristics and thus bring shear-extension coupling into the analysis. When 
laminates are unsymmetrically stacked, bending-stretching coupling is added to 
the analysis. Moreover, differences in elastic properties between fiber filaments 
and matrix materials lead to inplane-shear coupling. The high ratio of mplane 
modulus to transverse shear modulus makes the classical lamination theory, which 
neglects the effect of out-of-plane strains, inadequate for the analysis of multilayered 
composites. In such a case, one should use a theory which includes transverse shear 
deformation, an example being the first order shear deformation theory (FSDT). 

Thus, we need a displacement field that would be able to capture the existence 
of the various coupling effects, which play a major roll in laminated composites. 
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The following displacement field for the first-order shear deformation beam theory, 
in the defined coordinate system of reference, can be used to take these couplings 
into account: 


U(x,y, z,t) = u{x,t) + y(3{x,t) + z<f>{x,t) - zg(x,y,t) (2.1a) 

V{x,y,z,t) = v(x,t) - zr(x ) t) ( 2 -lb) 

W{x,y,z,t) = w(x,t) + yr(x,i) ( 2 - lc ) 

where u{x , t) is the axial displacement, v(x,t) is the lateral displacement, w(x,t) is 
the transverse displacement, t) is the rotation of the transverse normals with 
respect to x, /3(x,t) is the in-plane rotation, t(x,£) is the twist angle, and g(x ) y,t) 
is the warping function to account for twist, bending, and extensional effects. All 
these displacements and rotations are measured at the midsurface. 

The generalized displacement vector is defined as follows: 

d T = | u(i,i) v(x,t) w(x,t) r(x,t) <j>(x,t) fi(x,t) ) (2.2) 


Further, we assume that it is possible to separate the displacement components 
of the laminated beam as products of time and spatial functions, assuming the same 
time function for each displacement, i.e., 


m 

m 

m 


U(x, y, z) = u{x) + y f3(x) + z <f>(x) - 2 g(x , y) (2.3a) 
V(x,y, z) = v(x) — zt(x) (2.3b) 

W(x,y,z) = w(x) + yr(x) 


(2.3c) 


2.1 CONTINUUM MECHANICS 


26 



Figure 2.2: Large deformation of a body from the initial configuration, C°, to the current 
configuration, C 1 

2.4.2 Displacement Gradients 

Figure 2.2 shows a general body in its initial configuration and in its current con- 
figuration (after deformation). Let the body in its undeformed configuration have 
a volume designated To, external surface area fi 0 , mass density po> and reference 
material points of the body to cartesian coordinates x 0 , yo, zo- Denote the deformed 
configuration with a volume Ti, external surface area Di, mass density pi, and ref- 
erence material points of the body to cartesian coordinates xi,yi,zi- Also, let the 
coordinate systems of the reference and current configuration coincide. 

The initial position of a point, P 0 , with coordinates {x 0 ,y 0 , to), is given by the 
position vector ro, and the current position of the same point, Pi, with coordinates 
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(xi,yi,zi), is given by the position vector 1 * 1 . The position vector ri is defined as 


ri = r 0 + R 


(2.4) 


where 



f > 

U{x,y,z,t) 


/ > 
Xo 


Xi 1 

R = i 

V{x,y,z,t) 

> r 0 = < 

yo 

> ri = < 

^ 1 


< W(x, y, z, t) j 


, 2 ° > 


2 i J 


(2.5) 


The derivatives of ri with respect to r 0 constitute the deformation gradient 
matrix, F, when arranged in Jacobian format: 


F = 


dx\ 

dxi 

dx i 


■ dU 

dU 

dU ~ 

dx 0 
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dz o 
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1 © 
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. dx o 

dyo 


(2.6) 


The determinant of the deformation gradient matrix is known as the Jacobian 
determinant and is defined as 

J = det[ F] (2.7) 

The displacement gradients with respect to the reference configuration are defined 
as follows: 


G = F — I = 


dx\ 


dx 0 
dy i 

dx 0 

dz i 
dx Q 


- 1 


axj 

dx i 


dU 
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1^) 

dyi x 
dyo 

dy i 


ay 

ay 

ay 

a* 0 


arco 

dy o 

az 0 

dz\ 

a«i j 

az 0 


aw 

aw 

aw 

dyo 


. ax 0 
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(2.8) 
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For the Total Lagrangian description, it is convenient to arrange the displacement 
gradients in vector form as follows: 

g T = { 9i 92 93 94 93 96 97 9s 9s } ( 2 - 9 ) 


Since the strains and displacements are referred in the reference configuration, 
we will take r = r 0 . Now, the displacement gradients for the assumed displacement 
field in Eq. (2.3) are: 


9 1 
9 2 
9 3 
94 


9 5 


9 6 


9 7 
9 8 
9 9 


dU(x,y, 

z,t) 

dx 


dV(x,y , 

z,t) 

dx 


dW(x,y, 

z,t) 

dx 


dU(x,y, 

z,t) 

dy 


dV(x, y, 

z,t) 

dy 


dW(x, y. 

,Z,t) 

dy 


dU(x,i /, 

z,t) 

dz 


dV(x,y, 

z,t) 

dz 


dW(x,y 

,z,t) 


dz 


du d/3 d<j> 

dx dx dx 

dv dr 

dx Z dx 
dw dr 

dx ^ dx 
= /? 

dy 

= 0 
= T 


= <f>-g{x,v) 


= —T 


= 0 


dg(x,y) 
Z dx 


(2.10a) 

(2.10b) 

(2.10c) 

(2.10d) 

(2.10e) 

(2.10f) 

(2.10g) 

(2.10h) 
(2. lOi) 


2.4.3 Strain Measures 

The strains associated with the displacement field defined in Eq. (2.3) are computed 
using the Green-Lagrange strains. These strains can be expressed in terms of the 
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displacement gradients, in the quadratic form, as follows: 


e t = hfg + -g T H<g 

(2.11) 

where the 9x1 vectors h, s and 9 x 9 matrices Hj’s are given in 

Appendix A. 

Here we use the Von Karman nonlinear strains, which assume small strains and 
moderate rotations, and these can be expressed in terms of the Green-Langrange 

strains as follows: 


1 2 

Cl — c xx g i 2 ^ 

(2.12a) 

1 2 

e 2 = e yy — g$ + - g 6 

(2.12b) 

e 3 = e zz = g$ 

(2.12c) 

e4 = 2 e yz = g 6 + 98 

(2.12d) 

e5 = 2e xz = fife + gi 

(2.12e) 

ee = 2 e xy = 32 + g* + 9 z 9 e 

(2.12f) 

Note that e 3 = e 4 = 0 for our displacement field. The strains in Eq. ( 2 . 12 ) can be 
expressed in terms of the midplane strains and curvatures as follows: 

“ ~ —.0 1 y, ^ - 0 

&XX E xx "T Z ™xx 

( 2 . 13 a) 

0 J 0 

e yy ~ e yy z K yy 

( 2 . 13 b) 

2 e xz = 7 ° xz + z K xz 

( 2 . 13 c) 

Z ^ xy 

( 2 . 13 d) 

where (e° x , e° yy , -f° z , l° xy ) are the membrane strains and « x , 

*° yy , < z , «x„) the 

bending strains (curvatures). Now, we express the midplane strain and curvature 
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terms in vectorial form, and further separate them into linear and nonlinear strains. 

(2.14) 


e° = £l + £nl 


where £l is the vector with midplane linear strains, 

— f L ^ L k L k L 1 ^ 

“ \ ^xxi ^yy’ 7xy > "'xx’ ^yy ’ ™xy > 'xzj 

and is the vector with midplane nonlinear strains, 

- - I e NL r NL k nl k nl y NL \ T 

^NL — \^xx i^yy » 7xy > ^yy ’ ^xy i 'xz / 

The nonzero midplane linear strains are 


II 

^ H 
Uj 

du dp 

Tx + y di 

(2.15a) 

II 

d<p dg{x,y) 
dx dx 

(2.15b) 

II 

V 8 

dw dr . 

~dx + ^ + V 'di ~ 9 ^ x ' y > 

(2.15c) 

II 

dv 

Tx + P 

(2.15d) 

K *v = 

dr dg(x, y) 
dx dy 

(2.15e) 


The nonzero midplane nonlinear strains are 


7: 




'yy 

7VL 

xy 


1 f/3w\ 2 , „ dm dr ( ft-Vl 

= + 2 !, fe 8 i + l ! 'toJ | 


1 ( \2 

= 2 (r) 

/ ckc dr 

- & + !, & |T 


(2.16a) 

(2.16b) 

(2.16c) 
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2.4.4 Warping Function 


Although many researchers have suggested various methods of approximating the 
warping distribution, here we use the simplest possible warping function. The func- 
tion is picked such that it satisfies the classical laminated plate theory (CLPT), 


which assumes that the transverse normals rotate such that they remain perpen- 
dicular to the midsurface after deformation. In the CLPT 7 = ^ an d this is 


satisfied by 


<p + dw /dx = 0 and 


V ^x ~ 9 ^ x,y ^ = 0 


Further adding a warping constant a, the warping function used here becomes: 


, . dr 

g(x,y) = ay^ 


(2.17) 


where a takes values of 0 or 1: when a = 0, warping is not included, when a 1, 
warping is considered. 

Thus using this warping function, the transverse shear strain and curvatures 
k zx and K zy , for the FSDT, are expressed as follows: 


7x2 = ^ + <t> + v( !-«)£: 


dw 
dx 

T d<p d 2 r 
= dx~ ay d ? 


dr 

dx 


K, 


xy 
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2.4.5 Stress Measures 


The stresses corresponding to the Green-Lagrange strains are the second Piola- 
Kirchhoff stresses (PK2). The three dimensional tensor in cartesian coordinates is 


Sxx 

Q 

u xy 

s xz ‘ 


S X x 

S xy 

Sxz 

Q 

Llyx 

Syy 

Syz 

= 

S X y 

q 

u yy 

Syz 

m S zx 

S Z y 

s zz . 


S X z 

Syz 

s zz . 


It can be shown that the PK2 stresses are linearly related to the Cauchy stresses 
as follows (Truesdell and Noll, 1965): 


S = S° + JF' 1 <tF _t 


(2.19) 


where S° are the prestresses, J the Jacobian determinant, F the deformation gra- 
dient matrix, S the PK2 stresses, and <r the Cauchy (true) stresses defined as 


cr = 


0 X x & xy O xz 

Oxy Cyy O' y Z 
O xz ®yz ® zz 


( 2 . 20 ) 


Not only must the transformations that are the motions for each of the mass- 
elements in an ensemble obey mass conservation, but the total mass of the entire 
body must be conserved (McDonald, 1996). Thus, 


P\dT\ = p 0 dT 0 


pi det{ F] dl o = Po dU o 


J = det[ F] 


<£T j 

dT o 


Po 
P i 
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where dr! and dT 0 are the volumes in the current configuration and reference 
configuration, respectively; p\ and po are the mass densities in the current and 
reference configuration, respectively. 

Assuming that isochoric deformation takes place (volume-preserving deforma- 
tion), J = 1. Also, we assume that the stresses in the reference configuration, 
S°, are zero. Lastly, recall that we restrict our analysis to large deformations but 
small strains. Under these assumptions, it can be shown that the PK2 and Cauchy 
stresses approach each other (Truesdell and Noll, 1965). Thus Eq. (2.19) reduces 
to 

S«<r (2.21) 

2.5 Laminate Constitutive Relations 


The relation between stress and strain is given by a constitutive equation. In 
this section we present the constitutive equations for linear hyperelastic laminated 
composite plates. 

2.5.1 Linear Hyperelastic Material Law 

Although the strain e zz was found to be zero, the present analysis will assume a state 
of plane stress (S„=0) and condense e« from the stress-strain relationship. The 
reduced material coefficients are expressed as Qij. Since in laminated composites 
each ply may have different orientation, the stresses are expressed in terms of an 
arbitrary angle 9 and the transformed plane stress-reduced elastic coefficients are 
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expressed in the 3 x 3 matrix Q (see Appendix B.l for details): 


Q 


Q 11 Q 12 Ql6 
Q 12 Q22 Q26 
Q 16 Q26 Q 66 


( 2 . 22 ) 


The PK2 stresses are expressed in the 3 x 1 vector S and the inplane strains are 
expressed in the 3x1 vector e: 



(2.23) 


Further assuming that the material has monoclinic symmetry with respect to the 
reference plane of the beam, the transverse shear strains are uncoupled from the 
inplane strains in Hooke’s Law. Thus the transformed stress-strain relationship 
takes the following form: 


S = Qe 


(2.24) 


The transverse shear stresses and shear strains are related by Hooke’s Law as fol- 
lows: 



Qu Q45 f 2e yz 1 

Q 45 Qs5 . i ‘4e xz j 


(2.25) 


These stresses when evaluated at the k th lamina are expressed with a superscript 
k. 
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2.5.2 A, B, D Matrices 


The integration of the stresses throughout the thickness leads to generalized stresses 
(force and moment resultants): 
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Qx = / S X z dz 


<0 

II 

/ Syx 

dz 

J -h/2 



J —h/2 



*5 1 xy dz 


S xy z dz 


where h is the thickness of the beam. Although the strains are continuous through 
the thickness, stresses are not, due to change in material coefficients through the 
thickness. Hence, the integration of the above stresses through the laminate thick- 
ness requires a laminawise integration (Reddy, 1997). 

The integration leads to the well-known matrices: the extensional matrix A, 
the extensional-bending coupling matrix B, and the bending stiffness matrix D. 
These are defined as 


/ hi 2 -Mam 

Qijdz ='^2Q ij (z k+ i - z k ) i,j = 1,2,4, 5,6 (2.26) 

V 2 fc=l 

/ h/2 / 2 ~2 \ 

Q v zdz U = (2.27) 

■h/2 fc=i V Z / 

rh/2 / 3 _ ~3\ 

l ’ j = 1 ' 2 ’ e (2 ' 28) 


Nla 


where 7V lam is the total number of plies considered. When considering symmetrically 
laminated composites, B is identically zero and the coupling between bending and 
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stretching vanish. However, for unsymmetrically laminated beams, the coupling 
cannot be ignored and it must be included in the analysis. 

2.5.3 Constitutive Equations for Transverse Shear Resul- 
tants 

The integration of the transverse shear stresses, Eq. (2.25), through the thickness 
of a laminated plate yields the following relation: 

f Qy 1 ^44 fei fcj ^45 fTyzl (2 29) 

j ^1 ^-2 -^45 ^2 ^55 _ \Tx2 j 

where k\ and k\ are the plate shear correction factors (Whitney, 1972). For the 
case of laminated beams Q y = 0, thus we condense out the shear strain 7 yz from 
Eq. (2.29). The constitutive relation for the transverse shear resultant simplifies to 

Qx = k\ ^55 — Txz (2.30) 

Let us define the equivalent bending stiffness due to shear as 

D„ s = (231) 

and redefine the shear correction factor as K a = k%. Thus the constitutive relation 
for the transverse shear resultant of the laminated beams considered here becomes 

Q x = K a D C55 7x2 ( 2 ' 32 ) 

Here the values of K a are assumed as 5/6, a value often used in the literature. 
However, a more rigorous treatment is given by Cohen (1978). This is discussed in 
more detail in a later section. 
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By integrating Eqs. (2.24) through the thickness, the basic constitutive relation 
can be expressed in terms of the A, B, D matrices. Further, by adding an addi- 
tional stiffness coefficient corresponding to the transverse shear modulus, by using 
Eq. (2.32), we get the following constitutive relations: 
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(2.33) 


In the analysis of one-dimensional structures, 

Myy = N yy = 0 


(2.34) 


Thus, the midplane strain e° yy and the bending curvature K° yy are condensed from the 
constitutive equations. This is done by rearranging and partitioning the bending- 
stiffness matrix: _ 


D c = 


Di,i Dili 
_Dnj Dii,h 

[0] 


[0] 


[DcJ 


(2.35) 


where the first submatrix in the bending-stiffness matrix is reduced as follows: 


Dr = Di i - Di,n D n j! Dh i 


(2.36) 
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This leads to the following equivalent bending-stiffness matrix for an unsymmetri- 
cally laminated beam: 


An 

Ai 2 

A13 

A„ 

0 

A 12 

A 2 2 

“A 2 3 

A 2 4 

0 

Ais 

A 2 3 

A 3 3 

A34 
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0 

K.D, 


(2.37) 


The details of the derivation and the coefficients of D c are given in Appendix B.2. 


Thus the laminate constitutive equations in Eq. (2.33) reduce to 
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(2.38) 


2.6 Nondimensionalization 


In many problems, we are interested in comparing the dimensionless response rather 
than the actual values. The nondimensionalization is of great help when we are 
comparing results with different properties. Two ways exist to nondimensionalize: 
(i) solve the problem dimensionally and nondimensionalize the response, or (ii) 
nondimensionalize the equation, resulting in a dimensionless response. In many 
problems, the latter approach is a more elegant one because it helps us perform 
parametric studies. Since in the literature very little attention has been given 
regarding the nondimensionalization of the laminated equations of motion, here we 
present a systematic way of doing so. 
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For the shake of convenience, let us define the following non-dimensionalized 


parameters: 


eh e 

r ' = h r2 = i nr2 = i 


Moreover, the coordinates x, y, and z are nondimensionalized as follows: 


(2.39) 



(2.40) 


2.6.1 Dimensionless Displacement Field 


The displacement field in Eq. (2.3) is expressed in terms of the midplane displace- 
ments and rotation. Thus we begin with the nondimensionalization of these dis- 
placements and rotations. Note that the twist rotation, r, the rotation of the 
transverse normals with respect to x , 0. and the in-plane rotation, 0, are all di- 
mensionless quantities because their values are given in radians. The axial, lateral, 
and transverse displacements are nondimensionalized as follows: 


u v 



w 



(2.41) 


Let us divide Eq. (2.3) by e and note that 

y _ y b _ y 
e be n r 2 


Thus the dimensionless displacement field for the first-order shear deformation 
beam theory becomes 


U{x,y, z) = u(x) + + j(j){x) 

z 

V(x,y,z) = v(x)-jt(x) 

W(x,y,z) = w{x) + ^-t(x) 


y z^ dr{x) 
nr 2 e dx 


(2.42a) 

(2.42b) 

(2.42c) 
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2.6.2 Dimensionless Strains 

The displacement gradients given in Eq. (2.10) are expressed in terms of the nondi- 
mensional quantities as follows: 


01 = 

du _ au 

dx dx 

(2.43a) 

02 = 

dV _ dV 
dx dx 

(2.43b) 

03 = 

'ik ** 
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II 
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(2.43c) 

04 = 

du dU 
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(2.43d) 
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dV dV 

dy dy 

(2.43e) 

06 = 

dW dW 
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(2.43f) 
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dU _ dU 
dz dz 

(2.43s) 

08 - 

dV _ dV 
dz dz 

(2.43h) 

09 = 

dW _ dW 
dz dz 

(2.43i) 


Thus dimensionless nonzero midplane linear strains can be expressed in terms of 
the dimensionless quantities as follows: 
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The dimensionless nonzero midplane nonlinear strains can be expressed in terms of 
the dimensionless quantities as follows: 


,NL 

XX 


j NL 

Ixy 


1 f / dw\ 2 2 y dui dr ^ f y 3t\ 2 1 

2 1 V dx ) ri r 2 dx dx \ r i r 2 dx ) J 

/ dw y dr\ 

\ dx T\ r 2 dx ) 


(2.45a) 

(2.45b) 


These dimensional strain quantities are related to the nondimensional strains as 
follows: 


£ L = 


£ L 



£ NL = £ NL 


,NL 


= 7 NL 


(2.46) 


2.6.3 Dimensionless Constitutive Equations 

The integration of the stresses through the thickness leads to the laminated consti- 
tutive equations, which is expressed in terms of the 3x3 extensional, extensional- 
bending coupling, and bending stiffness matrices as follows: 

N = A e° + B k° (2.47a) 

M = B£° + Dk° (2.47b) 

We nondimensionalize the stress resultants as follows: 

N ^ M 

Eyy h Eyy k ^ 


(2.48) 


2. 6. NONDIMENSION A LIZ A TION 


42 


Thus Eq. (2.47) can be written as 

N = A _ c B *° 

Eyy k Eyy k Eyy k l 

M _ B _q D k° 

Eyy ^ ~ Eyyh 2 * Eyy ht l 


(2.49a) 

(2.49b) 


Note that the underlined terms can be rearranged and expressed in terms of di- 
mensionless quantities, given by Eq. (2.39), as follows: 


B k° B h R° _ B J_-o 

Eyyk T ~ Eyy If £ ^ Eyy T X 

D k° D h R° _ D _1_ _o 

Eyy h 2 T ~ Eyy h 3 £ _ Eyy ^ V, ^ 


(2.50) 


Let us define the dimensionless 3x3 extensional, extensional-bending coupling, and 
bending stiffness matrices as follows: 


A _ A s _ 1 B 

Eyyk Vi Eyyk 2 


- 1 D 

" fi Eyy h* 


(2.51) 


In order to express the laminated constitutive equations in terms of the dimen- 
sionless extensional, extensional-bending coupling, and bending stiffness matrices, 
we multiply Eq. (2.49b) by l/rq. Then the dimensionless laminated constitutive 
equations are expressed as follows: 


N = Ae° + BK° (2.52a) 

-M = B e° + D k° (2.52b) 

r i 


The condensation of £° yy and R° yy can be done in the same way as explained in 
Appendix B.2. Similarly, the dimensionless constitutive relation for the transverse 
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shear resultant of the laminated beams here considered becomes 

Qx = = Kt Ass %x (2 53 ) 

tjyy tl 

Thus the dimensionless laminate constitutive equations in Eq. (2.38) becomes 
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(2.54) 


In this work, the dimensionless generalized stress vector will be defined as 
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and the corresponding dimensionless generalized strain vector as 
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(2.55) 


(2.56) 


where the above dimensionless linear and nonlinear midplane strains are given by 
Eqs. (2.44) and (2.45). 
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2.7 Transverse Shear Correction Factor 

The first-order Mindlin theory assumes that the transverse shear strain and thus 
the transverse shear stress are constant along the thickness direction. However, the 
actual variation of the transverse shear stresses is not constant and this can lead 
to significant discrepancies for unsymmetrically laminated beams (Madabhusi and 
Davalos, 1996). To overcome this problem, one can use the transverse shear strain 
energy of the first order theory and compare it to the transverse shear strain energy 
of a more accurate shear distribution. This comparison will yield an effective shear 
correction factor. 

A common practice for obtaining a more accurate transverse shear correction 
factor is to obtain the transverse shear stresses by directly integrating the equilib- 
rium equations with respect to the thickness coordinate and assuming the inplane 
stresses vary linearly in z. In doing so, constants are introduced which are nor- 
mally determined from the zero shear traction condition at the bottom surface of 
the laminate. Prom equilibrium point of view, the transverse shear stress should 
vanish at the top and bottom of the plate. Researchers such as Madabhusi and 
Davalos (1996) have derived a general expression for the shear correction factor of 
laminated rectangular beams with arbitrary lay-up configurations using this pro- 
cedure. However, their formulation does not guarantee the shear stress S xz to be 
zero at the top. 

As opposed to the above approach, Cohen (1978) provides a systematic method 
to compute the three shear stiffnesses A 44 , A\ 5 , and A 55 . Cohen’s formulation is 
based on Taylor series expansion about a generic point for stress resultants and 
couples which identically satisfy plate equilibrium equations. By applying Cas- 
tigliano’s theorem of least work, the statically admissible transverse shear stress 
distribution is determined. We use Cohen’s procedure, and apply the procedure to 
laminated beams, to calculate the shear correction factor for laminated beams. 
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2.7.1 Constant Transverse Shear Stress Distribution 


The transverse shear strain energy density, per unit area, for a constant transverse 
shear stress distribution is 


1 f h/2 

f^shear = n I ^xz 7 xz 

z J-h/2 


(2.57) 


Using Eqs. (2.32) and (2.25) we get 


U shear 


1 f k/ 2 

2 J-h/2 


Qx 


K. d 


dz = - 


Qx 


C55 


2 K s D Cm> 


rh/2 

J-h/2 


S r , dz = - 


1 Q 


2 K s D C5S 


(2.58) 


Although the shear strain 7 ° 2 is zero for our displacement field, we assume it as 
nonzero and further assume the shear stress S yz is negligible (S yz ~ 0). Thus the 
constitutive relation for the transverse shear, given in Eq. (2.25), reduces to 


S„ = 7« = 0 7 ;, (2.59) 

where G is the equivalent shear modulus and is piecewise constant in 2 (Reddy, 
1997). 


2.7.2 Actual Transverse Shear Stress Distribution 


The transverse shear strain energy density, per unit area, for the actual variation 
of the transverse shear stress through the thickness is 

1 r h / 2 1 f h/ 2 S 2 

f/ 8hear = \ S xz G~ l S xz dz=\ -f dz (2.60) 

2 J-h/2 1 J-h/2 G 

Now we proceed to highlight the steps suggested by Cohen (1978) to calculate 
the above integral. The transverse shear stress S xz is obtained by integrating the 
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three-dimensional equilibrium equation (Whitney, 1987): 

dSxz 
dz 

S XZ 

Note that for the sake of convenience we have substituted z for Q in the integrand 
of Eq. (2.62). Now we proceed to calculate the stresses S xx and S xy . To do so, we 
assume that the two-dimensional constitutive equations for stretching and bending 
can be expressed as the usual linear distribution of strain: 

e = e° + C «° ( 2 - 63 ) 


dS xx _ 
dx 

d r 


9S x y 

dy 


~ / S xx dC 


dx 


h/2 


d_ r 

dy J -h /2 


S x y dC, 


(2.61) 

(2.62) 


Now Eq. (2.24) becomes 

S = Q £° + C Q 

The integration of the above equation through the thickness leads to 

N = A e° + B k° 

M = Be° + D/c° 


(2.64) 


(2.65a) 

(2.65b) 


where A, B, D are the 3 x 3 extensional, extensional-bending coupling, and bending 
stiffness matrices. The inverted relation of Eq. (2.65) is given as 

e° = A* N + B* M (2.66a) 

k ° = B* N + D* M (2.66b) 


Substituting Eq. (2.66) into Eq. (2.64), we get 

= [a + cb] n+ [: 

f(i) 



s 


(2.67) 
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where S, N, M are 3x1 vectors, and f^, axe 3x3 matrices. The 3 x 3 A, B, 
D matrices are calculated as follows: 

D* = [D — BA -1 B] _1 
B* = -A -1 BD* 

A* = [I - B* B] A -1 

A = QA* B = QB* D = QD* 

Now we substitute Eq. (2.67) into Eq. (2.62). Note that only the matrices ^ and 
f( 2 ) are functions of £. Thus these can be integrated in the thickness direction and 
expressed as 

F (1) = f Z f ( 1 ) (C)d£ F (2) = [ f ( 2 ) (C)^C 

J-h/2 J-h/2 

As a consequence, Eq. (2.62) becomes 

S xz = (ViV N X X + F\2 Nyy + E\[i N X yj 

- (Fj (2) M xx + Fg> Myy + Fg } M X y) 

and taking the partial derivatives, we get 

S XZ = - (Fl? N XX , X + F$ Nyy, X + Fg> N X y, X ) 

- (F^ M XX ,y + Fg } Myy,y + F^ M X y ^ 


( 2 . 68 ) 
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The partial derivatives of the stress resultants are calculated by satisfying the equa- 
tions of equilibrium (Whitney, 1987): 


NxXyX xy,y 0 

(2.69a) 

^Xy,X H" Nyy y y = 0 

(2.69b) 

Qx,x Qy,y 0 

(2.69c) 

Mxx,j “f" M X y,y Qx 

(2.69d) 

M X y } X “t" Myy^y Q y 

(2.69e) 

The above equations are satisfied by expanding the stress measures in Taylor series 

about the reference axis and only keeping the first-order terms in 
N and M: 

the expansion for 

N X x ci ^ ~b dj y 

(2.70a) 

N yy = ai x + 6i y 

(2.70b) 

N xy = -fei X - Cl y 

(2.70c) 

M xx — ( c 2 + Qx) x + d 2 y 

(2.70d) 

M yy = a 2 x + (b 2 + Q y ) y 

(2.70e) 

M xy = -b 2 x-c 2 y 

(2.70f) 

For the case of laminated beams we assumed that N yy = M yy ■ 
this assumption it can be shown that b\ = b 2 = 0. 

= Q y = 0. Under 


Substituting Eq. (2.70) into Eq. (2.68), and using the assumptions of laminated 
beams, we obtain a statically admissible transverse shear stress distribution. The 
result can be expressed as follows: 


— S xz = Qi £i + C*2 £2 + P Qx 


(2.71) 
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where & is the 2 x 1 coefficient vector defined as 

$7 = {ci di } i = 1,2 

The 1 x 2 vectors Qi,a 2 are functions of z and contain the polynomials correspond- 
ing to the coefficients &. The lxl vector (3 is also a function of 2 . These are 
obtained using the symbolic processor MATHEMATIC A 1 and are given as 


«i = 


h An - h 2 Bn z Bn 
- t- 2 A n g 2 


T 

<*2 = 


^ h 2 D33 2 2 D33 h B33 

8 2 2 ~ 

h 2 Dj 3 2 2 D 13 fiBia 


2B 


33 


8 


+ 


+ 


+ 2 B 


13 


/3 = 


/i 2 Dn . 2 2 Du . fiBn 


+ 


- + 


+ 2 B 


11 


8 2 ' 2 

The substitution of Eq. (2.71) into Eq. (2.60) gives the transverse shear strain 
energy per unit area. Cohen (1978) suggests expressing this transverse shear strain 
energy density as 


f/shear = ^ X T A X + X T 1 Q + ^ Q T C Q 


(2.72) 


where A is a 4 x 4 matrix defined as 


An A12 
A -12 A -22 J 


where A 


V2 

[ a i Q j 

J G 

- h / 2 


dz 


1 A registered trademark of Wolfram Research, Inc. 
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B is a 4 x 1 vector defined as 

f Bi 

B = < _ 

1 b 2 

C is a 1 x 1 matrix defined as 


h / 2 

a . 1 3 

where B; = / *„ dz 


V '~ I 

-N 2 


G 


hi 2 

c- f 


-I 

-hjl 


G 


is a 4 x 1 vector defined as 


rT = U? ^2 } = {ci di c 2 d 2 } 


and Q is a 1 x 1 vector defined as 


Q = Q X 


Note that the integration for Ay, Bj, and C must be done laminawise. The details 
of the derivation and the coefficients of the above matrices and vectors are given in 
Appendix B.3. 

Since the shear stress distribution S xz is statically correct for all values of X, it 
follows from Castigliano’s theorem of least work (Timoshenko and Goodier, 1970) 
that the proper values of X are those which minimize U e hear- The minimization of 
Eq. (2.72) with respect to X are the stationary values of the transverse shear strain 
energy density. In other words, 

— hear = AX + BQ = 0 => AX = — ® Q 


Two possibilities exist: when A is nonsingular and when A is singular. Cohen 


2. 7. TRANSVERSE SHEAR CORRECTION FACTOR 


51 


(1978) provides the complementary transverse shear strain energy for both of these 
cases. 


(i) A is a nonsingular matrix: |A| ^ 0 

Rehear = ^ Q T [C - B T A’ 1 B] Q (2.73) 

(ii) A is a singular matrix: |A| = 0 

First calculate the eigenvalues and eigenvectors of A. Let $ be the orthogonal 
4x4 matrix whose columns are orthonormal eigenvectors of A. Thus, 

$ t A* = A (2.74) 

where A is the 4x4 diagonal matrix containing the eigenvalues and the f th diagonal 
element corresponds to the ? th column of Now let us suppose there are p zero 
eigenvalues for p < 4. Then we eliminate the columns of the eigenvector matrix 4* 
corresponding to these zero eigenvalues. This will lead to a 4 x (4 — p) matrix $ p . 
Now, we obtain the new (4 — p) x (4 — p) diagonal matrix A p by retaining only the 
(4 — p) nonzero diagonals. 

Next we proceed to calculate A -1 (a 4 x 4 matrix), known as the natural inverse 
of A, as follows: 

A -1 = 3> p A“ 1 <l>p (2.75) 

Then the solution proceeds as before: 

U s hear = [C - ® T A _1 ®] Q 


(2.76) 
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2.7.3 Energy Equivalence Principle 

Note that in general the transverse shear strain energy density due to the actual 
transverse shear stress distribution can be written as 

(C-B t A-‘B) (2.77) 

Now, equating Eq. (2.77) and Eq. (2.58), we get the following: 

K • = D„ (C-B'A-'B) <2 78) 


2.8 Generalized Principle of Virtual Work 


The equations of motion for the deformed body in Fig. 2.2 can be developed using 
the principle of virtual work. Virtual work is defined as the work done by actual 
forces in displacing the body through virtual displacements that are consistent with 
the geometric constraints imposed on the body and possess sufficient continuity 
to compute virtual strains. The principle of virtual work states that a body is 
in equilibrium if and only if the virtual work of all forces is zero for any virtual 
displacement. Moreover, the principle of virtual work can be extended for dynamic 
analysis and nonconservative loads. Thus in this work, we use the generalized 
principle of virtual work (GPVW), or the extended Hamiltonian’s principle, and 
define it as follows: 



+ ^VVi n t 


- <WC - 6W n C } dt = 0 


(2.79) 


where £W ex t is the virtual work done by external forces, <5Wi n t the virtual work done 
by internal force, 5fC the virtual kinetic energy, and 5W nc the virtual work done 
by nonconservative forces. In this section, we derive this expression and further 
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nondimensionalize them. 


2.8.1 Virtual Kinetic Energy 


The virtual kinetic energy is given by 

8IC = JJJ p [u8U + V6V + W5w}dT (2.80) 

r 

Substituting the displacement field described in Eq. (2.3), we get 

6/C = Iff p\ (u + y + z<j> — zyf'^j ^6u + y8f3 + z84> — zySf'j 

r l ( 2 . 81 ) 

+ (v — z t) (Sv — z 8f) -T (w + y t) ( 8w + y 6f ) ^ dT 


Now, expanding the above equation, we get: 


6/C 


= JJJ pj (u + y/3 + z<i>- zySr'^) 8u 

+ (u + y/3 + z<p-zy 6f y 6/3 
+ (u + y/3 + z<j) — zy Sf'^j z8(f> + (v — zi) 8v 
+ (u + yfi + z<fi- zyT ,S j (—zySf 1 ) 

+ (v — zt) (8v — z 8f) + (w + yf) ( 8w + y 6f ) 

+ ( w + yi ) 8w+ (—zv + z 2 f + yw + y 2 f) 6f| dY 


The underlined terms are ignored because it is assumed that the contribution of the 
warping function to the overall kinetic energy is very small. Also, recalling that, 
in the GPVW, the virtual kinetic energy is integrated over the time domain, we 
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integrate by parts with respect to time to express the virtual velocities in terms of 
virtual displacements. Thus, 


f 


T l h/2 6/2 


SKdt 


1 c uj * vf * , 

— — / / p < (il + z(pj 8u+ (v - zt) 8v + wS 

0 0 -h/2 — 6/2 ' 

+ y 2 (ii — zf_ j 6/3 + y {t 5w + id <5r} 

+ (z u + z 2 <Pj 5(j)+ (-z v + y 2 f + z 2 f) St 

+ y | (j3 - z f 6u+ (u + z(fj 5(3 + (j3 — z f Sep} j dy dz dx 


dt 


+ /M 


l t=T 


t= 0 


where p is the mass density, h the total thickness, b the total width, and t the 
total length. Note that the underlined terms in the above expression vanish as they 
are integrated through the width. The double underlined term is ignored because 
its contribution is assumed negligible, i.e., zr f '{3. The boundary term f(t) 
is evaluated at two different times and vanishes because the Hamilton’s principle 
assumes that the actual dynamic path coincides with the varied path at the two 
time instants t = 0,T. The above equation can also be written as 


T T i 

5tCdt = - 
0 0 0 


1 

I 



I 0 u + h<pj 6u + ( I 0 V - Ilf) 8v + (low) 5w 


(2.82) 


+ 


(— I\v + J x t) St + (^hu + 5<f>+(^Jy'(3j <5/3>d:r<if 


where Io, I\, h are the mass moments of inertia, J x is the polar mass moment of 
inertia, and J y is the inplane mass moment of inertia. The integration of the above 
mass inertias through the laminate thickness requires a laminawise integration. 
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These are defined as follows: 


rh/2 rb/2 

I 0 = / pdydz 

J-h/2 J-b/2 

= by~]p k {z k + 1 - Zk) 

k = 1 

(2.83) 

rh/2 rb/2 

h= / pzdydz 

J-h/2 J-b/2 

Mam / 2 __ ^2 \ 

1 V 7 

(2.84) 

rh/2 rb/2 

h= / pz 2 dydz 

J-h/2 J-b/2 

Mam / 3 r 3 \ 

fc=l x 7 

(2.85) 

rh/2 rb/2 

J x = p (y 2 + z 2 ) dy dz 

J-h/2 J-b/2 

b 2 T _ 
= 12 /o + /2 

(2.86) 

rh/2 rb/2 

J y = / py 2 dydz 

J-h/2 J-b/2 

_ b 2 
"l2 /o 

(2.87) 


where JV, am is the number of plies. Kapania and Raciti (1989a) assumed the polar 
mass moment of inertia as J x — b 2 io/12. However, for rectangular cross sections 
this is not true. Thus, here we include the complete expression for J x . 

Since the virtual kinetic energy in Eq. (2.82) is only a function of acceleration, 
let us define the contribution of the virtual kinetic energy in Eq. (2.79) in terms of 
the virtual work done by inertial forces as follows: 

f T SKdt= f SW ineT dt = -f [ Sd r Mddxdt ( 2 . 88 ) 

Jo Jo Jo Jo 

where <5W iner is the virtual work done by inertial forces. Thus Eq. (2.79) becomes 

[ T {-<SW e xt + 6W mt ~ (SVViner - 5W nc } dt = 0 (2.89) 

Jo 

Now we proceed to nondimensionalize the virtual kinetic energy. The mass 
moments of inertia, the polar mass moment of inertia, and the inplane mass moment 
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of inertia axe nondimensionalized as follows: 


l o = 


Ii o 
pbh 


h = 


h 


pbh z 


h = 


pbh* 


(2.90) 




1 T 


Jx ~pb 3 h 12 


= — 7 0 + r 2 / 2 Jy — 


Jy _ 1 r 

— — in 


pb 3 h 12 


The limits of integration are changed as follows: 


(2.91) 


x = £x 


dx = £dx 


Thus Eq. (2.82) is expressed as follows: 


<5W lt 


= — j i^(pbh£ Iqu + pbh 2 Ii<j>) £5u 


+ (pbht I 0 v - pbh 2 hr) £5v + (pbh£ I 0 w) £5w 
+ (-pbh 2 £I\V + pb 3 hJ x r) 5 t + (pbh 2 £hu + pbh 3 1 2 (p) <50 

-I- (pb 3 h J y 0) £ dx 


Now we factor out pbh£ 3 to get 


<5Wj ner = — pbh £ 3 j | (j 0 u + j h + (l 0 v - jh Sv + (I 0 ui ) 6w 


+ [- t jhv+ Q) J x rj 6f + 




h<t> <50 


+ 11^1 J< 


» & \ <5/3 1 


dx 


Using the dimensionless quantities defined in Eq. (2.39), and further expressing the 
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result in matrix form, we get 


6W mer = -pb hP / M T Mil dx 


/■ 

Jo 


<5W ir 


(2.92) 


where 


6<P 


-{ 


5u 5v 6w 6 t 


M = 


lo 

0 

0 

0 

-h 

n 

o 


o 

/o 

0 

-h 0 

ri 


0 

0 


S(t> 8(3 } 


0 
0 
l o 


0 

0 


0 

--/l 

ri 

0 


2 2 
r l r 2 


0 

0 


-/l 

n 

o 

o 

0 


0 

0 

0 

0 

0 


2 „2 u v 

r l r 2 


(2.93) 


(2.94) 


It should be noted that the mass of the system is positive and that M is a symmetric 
and positive definite matrix. 


2.8.2 Virtual Work done by External Forces 


The virtual work done by external forces, in the absence of body forces, can be 


written as 


8W t 


ext 


= b e 



tj Sdj dQ, 


(2.95) 
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where tj are the dimensionless forces acting on the surface of the structure and 6dj 
the dimensionless virtual displacements. 


2.8.3 Virtual Work done by Internal Forces 

The virtual work done by internal forces can be written as 

SW M = jjj SjdejdV (2.96) 

r 

where Sj are the PK2 stresses and are energetically conjugate to e 3 , which are the 
Green-Lagrange strains. The virtual work done by internal forces is expressed as 


r t rb / 2 




int 


= [ [ {N xx 6e 0 xx + N xy 8 1 l y + M xx 5K° xx + M xy 5n° xy 

Jo J-b/2 [ ( 2 . 97 ) 

+ Qx h° x \dydx 


where I is the total length of the beam, and b the width of beam. 

Now we proceed to nondimensionalize the virtual work done by internal forces. 
The limits of integration are changed as follows: 


= lx 


dx = Idx y = by 


dy = bdy 


Thus Eq. (2.97) becomes 

rl rl/2 


sw \ 


int 


rl 1 _ _ l 

= J J lEyyhNxxS^ + EyyhNxyS^ly + Eyyh M XX “ <5 ^ 


1 


,1 


+ Eyy h 2 M X y ~ 8R° xy + Eyy h Q X 5^ ^ ^ ^ dX 
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Let us factor out E yy hb£ from the above equation: 


a: 


rl />l/2 ( _ - h _ 

6W mt = Eyy hb I I / { N xx 8t xx + N xy &f xy + M xx - Sk° xx 

-1/2 I 


+ Mxy J SK° xy + Q x Sj° xz dy dx 


} 


h 

and substitute — for r\ to get 


6W int = E yy hb£ 


1 r 1/2 


a 


M, 


1/2 


N xx 5? xx + N xy 6 ? xy + —6R: 


XX £ — o 
XX 


n 


M, 


+ ^5K° xy + Q x 6r xz }dydx 


n 


4 . 

xz ( 


( 2 . 98 ) 


Expressing the above equation in terms of the dimensionless generalized stress 
vector, Eq. (2.55), and dimensionless generalized virtual strains, Eq. (2.56), we 


get 


-o 

Pi 
r Pi 

II 

e 

< <5e T T i dy dx 

(2.99) 

*/ u %/ - 

-I/ 2 t J 



<*W int 



2.8.4 Virtual Work done by Nonconservative Forces 

The virtual work done by nonconservative forces is calculated using the definition 
of virtual work, 

<5W nc = 8d r Q (2.100) 

where Q is the vector of nonconservative forces and <5d the corresponding virtual 
displacements. 

The types of nonconservative forces considered here are shown in Fig. 2.3. Thus, 
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Figure 2.3: Laminated beam subject to a follower load 
for a follower load, Q becomes 

fb/2 _ 

Q t = i N xx cos 6 0 N xx sin# 0 0 0 \ dy (2.101) 

J-b / 2 L J 

<5d T = | 8u 6v 5w 8t 6(j) 5(3 j (2.102) 

In general, the tangential follower force remains parallel to the midsurface, i.e., 

e = <t>. 

Let the nonconservative axial load be P = N xx b. Then in order to nondimen- 
sionalize the virtual work done by nonconservative force, let us express Eq. (2.101) 
in terms of nondimensional quantities: 

£W nc = 5d T Q = N xx b cos <f> £ 5u + N xx b sin 4> i 5v 

Thus the virtual work done by nonconservative force is expressed as follows: 

<5VV nc = P t (cos 4> 5u + sin <f> 5v) 

' ^ 

5Wnc 


(2.103) 
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2.8.5 Equations of Motion 


Applying the GPVW, we get the following: 


J* { E yy h b i 6W int + pbh£ 3 5 W iner ~ P l 5W nc - b i <5W ext } dt = 0 

Now let us divide the above equation by E yy hb£ and multiply by (£/ h) 2 to get 
rT , pf ^ £ 2 


jfG 


SW mt + 4—, 5W i 


Eyy h 2 


iner , , 2 

Eyy hi l 


Eyy h 2 


5W ext }dt = 0 (2.104) 


2.9 Summary 


Since the noncognitive uncertainties are usually known in the structure’s reference 
configuration, the Total Lagrangian description is used because its reference con- 
figuration seldom changes. 

The assumed displacement field takes into account the various couplings that af- 
fect laminated composites: shear-extension coupling, bending-stretching coupling, 
and inplane-shear coupling. The present formulation also includes a warping func- 
tion. 

The Green- Lagrange strains are derived for the given displacement field. The 
stresses corresponding to the Green-Lagrange strains are the second Piola-Kirchhoff 
stresses. By assuming that isochoric deformation takes place, stresses in the refer- 
ence configuration are zero and, considering only small strains, the PK2 and Cauchy 
stresses coalesce. 

The laminated constitutive law is expressed in terms of the extensional matrix 
A, the extensional-bending coupling matrix B, and the bending stiffness matrix 


2.9. SUMMARY 


62 


D. Since a one-dimensional analysis is being considered, the load conditions are 
prescribed and the in-plane strain e yy and the bending curvature K yy are condensed 
from the constitutive equations. This leads to the constitutive equations used 
throughout this dissertation. 

Finally, the equations of motion for the Total Lagrangian description using the 
generalized principle of virtual work were derived. Further, these equations are 
expressed in terms of nondimensional quantities. 

In the next chapter, we proceed to discretize these equations. 


Chapter 3 

A Shear Deformable Laminated 

Beam Element 


Because aircraft structures consist of many irregular continuous components, 
an analytical solution is almost impossible for these structures. Moreover, as men- 
tioned before, composite materials have only increased the complexity of the solu- 
tion. The finite element analysis has been widely used in solving such problems. 
The finite element formulation can be obtained by discretizing the nondimensional 
form of the generalized principle of virtual work presented in chapter 2. 

In this chapter, we present the formulation of a twenty-one degree of freedom 
laminated beam element that takes into account the existence of various coupling 
effects, which play a major roll in laminated composite materials. This element is 
valid for the static and dynamic analysis of both symmetrically and unsymmetri- 
cally laminated beams. 


3.1 Discretized Continuum Mechanics 

In this section, we present the discretization of the beam using the Total Lagrangian 
formulation. 
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Figure 3.1: A twenty-one degree-of-freedom laminated beam element 

3.1.1 Displacement Field 


Since all quantities used here are nondimensional, we drop the overbar from the 
displacements in Eq. (2.42). Thus the dimensionless displacement field for the 
first-order shear deformation beam theory is rewritten as 

u(x,y,z) = (3 - la) 

V{x,y,z) = v{x)-jt(x) (3-lb) 

W(x,y,z) = w(x) + —t{x) (3.1c) 

n j *2 

where u is the midsurface axial displacement, v is the midsurface lateral displace- 
ment, w is the midsurface transverse displacement, (p is the midsurface rotation of 
the transverse normals with respect to x, (5 is the midsurface in-plane rotation, r 
is the midsurface twist angle, and a denotes the warping constant taking values of 
0 or 1. Here we define the ratios ri and as 



(3.2) 


where £ e , h, and b are the length, thickness, and width of the beam element, re- 
spectively. 
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In order to develop a shear deformable laminated finite element beam, one 
should avoid shear-locking (Reddy, 1993). Two methods exist to alleviate this 
problem: reduced integration (RI) and consistent interpolation elements (CIE). In 
this work, we use the CIE because this approach requires fewer degrees of freedom 
than the RI approach. Thus the orders of the interpolation functions are chosen to 
be consistent in the computation of the strains. 

To alleviate the shear-locking problem, the one-dimensional element is derived 
such that it consists of four equally spaced nodes and a node at the middle as shown 
in Fig. 3.1. The dimensionless nodal displacements measured at each node are: (i) 
at the exterior nodes ( nodes 1 and 2), axial displacement u, lateral deflection v, 
transverse deflection w, rotation of the transverse normals <fi, in-plane rotation /3, 
and twist angle r; (ii) at the two equally spaced interior nodes ( nodes 3 and 4), the 
torsion r, and the derivatives of v, w with respect to x\ (iii) at the middle node in 
the interior ( node 5), axial displacement u, and rotations (j> and (3. All these nodal 
displacements are measured at the midsurface and are expressed as follows 

q t = {ui,Vl,Wi,Ti,<f>i,Pl,U2,V2,W2,T2,<t>2,(32, 

u 3 , v 3 , v 4 , W 3 , w 4 , r 3 , r 4 , <t>3, fc} T 

The length and width of the element are nondimensionalized to unity, i.e., 

x 

0 < x < 1 where x = — , 

t e 

ii - y 

-- < y < 2 where y = 5’ 

and £ e , b are the actual length and width of the element, respectively. 

The midplane displacements v and w are chosen to obey a cubic polynomial of 
the form 

9t( f ) = y 

i = 0 


(3.4) 
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where s are coefficients to be determined using the following boundary conditions: 


<7i =Qt{x) 


x=0 


<72 =Qt{x) 


X=1 


qs = 


dqt 

dx 


x— 1/3 


q 4 = 


dqt 

dx 


i=2/3 


(3-5) 


where qt represents v or w. After substituting Eq. (3.4) into Eq. (3.5), the constants 
a, are obtained in terms of nodal displacements q,- Thus Eq. (3.4) can be rearranged 
and expressed in terms of the newly defined Lagrange polynomials as follows: 


q t (x) = (1 -4x + 9x 2 - 6x 3 ) qi + (4x-9x 2 + 6x 3 ) q 2 

y ^ ✓ v s / ^ 

N x {x) N 2 {x) 

+ (— 3x 2 + 3x 3 ) <73 + (3x - 6x 2 + 3x 3 ) q 4 

V ' N v 

N z (x) N a {x) 

The above shape functions are shown in Figure 3.2. The fact that cubic interpola- 
tion functions for w are used suggests that results will be close to those by an exact 
method. Moreover, it should be noted that by using these elements one ensures 
continuity in the displacements but not in their derivatives with respect to x. 

Although the CIE approach suggests to use a quadratic polynomial for the 
angle of twist, here we found that it does not produce good results when warping is 
included. One reason for this behavior could be that the inclusion of the warping 
function brings into the picture a second derivative in the curvature k xx and this 
requires a third-order polynomial to express the twist angle. Thus the twist angle t 
is chosen to obey a cubic behavior and it is represented using Lagrange interpolation 
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Figure 3.2: Nondimensional plot of the shape functions used to approximate the midplane 
displacements v and w. 

polynomials: 



$ 3 (x) M*) 


The above shape functions are shown in Figure 3.3. 


The axial displacement u and the shear rotations 4> and (3 are assumed to obey 
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Figure 3.3: Nondimensional plot of the shape functions used to approximate the midplane 
twist angle r. 


a quadratic polynomial of the form 

Tft(x) = (3.8) 

•=o 

where ry t represents & or 3. Using quadratic Lagrange interpolation polynomials, 
the following polynomials are obtained: 
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Figure 3.4: Nondimensional plot of the shape functions used to approximate the midplane 
axial displacement u and the shear rotations <j) and (3. 

The above shape functions are shown in Figure 3.4. 

Thus the deflection behavior of the beam element for the first-order theory is 
described as follows: 



u(x) 

= U\ + ^2 u 2 + ^3 ^3 

(3.10a) 

— 

v(x) 

= N\Vi + N 2 V 2 + N3 Vz + N 4 ^4 

(3.10b) 


w(x) 

= N\ W\ + U ) 2 + N% W 3 + N 4 W 4 

(3.10c) 

— 

t(x) 

= $lTi + <&2 7*2 + ^3 7*3 + $4 7*4 

(3.10d) 


4>{x) 

= 0i + \&2 02 + ^3 03 

(3. lOe) 

— 

0{x) 

= ^1 Pi + ^2 /?2 + ^3 /?3 

(3. lOf ) 
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The above equation can be rearranged in matrix form as follows: 


d = Nq t (3.11) 

where N is the dimensionless shape function matrix, and q t is the vector containing 
the dimensionless nodal displacements defined in Eq. (3.3). 


3.1.2 Strain-Displacement Relation 

Since the dimensionless virtual displacements are defined in the same space of 
functions as the finite element space of functions, the dimensionless virtual dis- 
placements corresponding to the dimensionless displacements defined in Eq. (3.11) 
are 

<5d = N<5q t (3-12) 

As a consequence, the total virtual generalized strains can be expressed in terms of 
the virtual linear and nonlinear dimensionless strains: 


5e = 8s l + 

(3.13a) 

c de c Del r , de N L x „ 

6e - — <5q t - _ 5q t + „ dq t 
«9q t oq t dq t 

(3.13b) 

= Bl <iqt + [BNL(qt)] 

(3.13c) 

= [Bl + B N L(qt)] 

(3.13d) 
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where Bl is the dimensionless strain-displacement matrix independent of the dis- 
placements and Bnl is the dimensionless nonlinear strain-displacement matrix lin- 
early dependent on the displacements. These matrices combined lead to the strain- 
displacement matrix, 

B s d = t . — = Bl + B NL (q t ) (3-14) 

dq t 

where qt is the vector containing the dimensionless nodal displacements defined in 
Eq. (3.3). 

3.2 Discretization of the Generalized PVW 


Generally, the finite element formulation is established in terms of a weak form 
of the partial differential equations under consideration. In solid mechanics, this 
implies the use of the principle of virtual work, or for dynamic and nonconservative 
problems, GPVW. Recall, from chapter 2, the equations of motion using GPVW 
were given by 

T + = ° (3 ' 15) 

The discretization of GPVW over the domain leads to the element tangent 
stiffness matrix, mass matrix, and force vector. 


3.2.1 Symbolic Computation 

The development of symbolic algebraic languages has made it possible to perform 
complex algebraic manipulations associated with various structural analysis prob- 
lems. The use of these symbolic manipulators reduces the possibility of making 
errors and considerably reduces the tediousness of solving complex problems by 
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hand. Moreover, the ability to develop analytical expressions for the tangent stiff- 
ness matrix avoids numerical integration errors and problems. 

In this work, the tangent stiffness matrix, the mass matrix, the internal force 
vector, and the external force vector are all obtained using the symbolic processor 
MATHEMATICA 1 . By obtaining these matrices and vectors analytically, the CPU 
time is reduced. This greatly saves CPU time for the Monte Carlo Simulation as 
one no longer has to perform numerical integration for each case. 

A very useful feature available in the symbolic processor MATHEMATICA is 
that it has the capability of providing the obtained analytical expressions in FOR- 
TRAN form. Thus, once the analytical expressions are obtained, these are easily 
implemented into FORTRAN 90. 

It would be even better to solve the entire problem using MATHEMATICA; 
however, it is time consuming and almost impractical for nonlinear problems. Com- 
puter programs such as FORTRAN are time-efficient in numerical calculation. 


3.2.2 Mass Matrix 


The virtual work done by inertial forces was defined in section 2.8.1 as 


5W ir 


= I «d T M3<a 

0 

Using the definition for d, Eq. (3.11), the above equation becomes 

= s q jym 


(3.16) 


SWu 


P MN dx 


(3.17) 


M c 


1 A registered trademark of Wolfram Research, Inc. 
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and M e is the dimensionless element mass matrix: 

M e = / N t MN dx (3.18) 

Jo 

where the length of the element is nondimensionalized to unity, N is the dimen- 
sionless shape function matrix, and matrix M is given by Eq. (2.94). 


3.2.3 External Force Vector 


The virtual work done by external forces, in the absence of body forces, can be 
written as 


<5W 


ext 



<5d T t dfi, 


(3.19) 


where tj are the forces acting on the surface of the structure and 5dj the dimen- 
sionless virtual displacements defined in Eq. (3.12). Using the definition for d, 
Eq. (3.11), the above equation becomes 

<5W ex t = JJ N T tdf} (3.20) 

h 

' 

F e 


Thus, the element external force vector is 

r 1 z‘i/2 

F e = / / N T t dydx (3.21) 

Jo J- 1/2 

where the width and the length of the beam element are nondimensionalized to 
unity. 
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3.2.4 Loading Stiffness Matrix 

The virtual work done by nonconservative forces is usually written in terms of the 
nodal displacements, i.e., 

<5W nc = &<ll fnc ( 3 - 22 ) 

where f nc are the dimensionless nonconservative force function of the nodal displace- 
ments defined in Eq. (3.3) and are defined for the global system but not locally, and 
Sctf are the dimensionless virtual nodal displacements corresponding to the nodal 
displacements corresponding to the force terms in the vector f nc . The dimensionless 
loading stiffness matrix is then obtained as follows (Argyris and Symeondis, 1981): 

K L = ^ (3.23) 

oq t 

The above matrix for nonconservative loading is unsymmetrical and is zero conser- 
vative loading. If the load is a distributed follower load, such as pressure, it can 
be included on an element basis. However, this matrix is obtained from the global 
nonconservative force vector. In other words, the loading stiffness matrix is added 
to the global tangent stiffness matrix (after assemblage) and not locally. 

Let us derive this matrix for the case of the follower load discussed in section 
2.8.4. We will start with Eq. (2.100): 

8W nc = 8d T Q (3.24) 

X— 1 

Note that the above expression is only evaluated for the element in contact with the 
follower load. Thus, the virtual displacements in the above expression are obtained 
by Eq. (3.12) as follows: 

<5d T =6cgN T (3.25) 

X=1 x=l 

where N T | , is the transpose of the shape function matrix evaluated at the second 
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node of the element, 



0 

0 

0 

0 

0 

0 

1 
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0 
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0 

0 

0 

0 

0 

1 

n| 

0 

0 

0 

0 

0 

0 

0 

0 

\x=l 
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0 

0 

0 

0 

0 

0 

0 

Nl 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 


0000000000000 

0000000000000 

1000000000000 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

(3.26) 

call the dimensionless 


0 1 0 0 0 0 0 0 

0 0 1 0 0 0 0 0 

0 0 0 1 0 0 0 0 

and 5qJ is given by Eq. (3.3). For sake of simplicity let us 

shape function matrix evaluated at the second node of the element as Nl. 


Now to calculate the dimensionless nonconservative force Q in terms of the gen- 
eralized nodal displacements, we use Figure 3.5. This figure shows the deformation 
of the last element at which the follower load is applied. The load rotates with 
the rotation at the tip of the cantilevered beam-column. Further we assume small 
rotation at the tip, 

cos 02 ~ 1 sin 02 ~ 02 (3.27) 

where 02 is the rotation of the tip. Also, note that the direction of 0 2 is taken 
opposite to the one assumed in out finite element formulation, as shown in Figures 
3.1 and 3.5. Thus, the dimensionless load Q is then expressed as follows: 


\ 

1 


0 ' 

0 


0 

0 

< 

> + < 

—02 

0 


0 

0 


0 

> 0 , 


0 


(3.28) 
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We can also express the above expression in matrix form: 


Q = Qi + Q 2 Qt 


(3.29) 


where 


q I = {«i V\ w l r l <Pl Pi u 2 v 2 w 2 r 2 <t>2 P2 u 3 v 3 v 4 w 3 w 4 r 3 r 4 4*3 @3^ 


Ql = { 


1 0 0 0 0 0 


} 


Q 2 = 


0000000000 

0000000000 

0000000000 

0000000000 

0000000000 

0000000000 


00000000000 

00000000000 

-10000000000 

00000000000 

00000000000 

00000000000 


Thus Eq. (3.24) becomes 

<5W nc = Sq[ N£ {Qi + Q 2 qt} (3.30) 

and the dimensionless elemental nonconservative force is then defined as 


^ c = nJQi + nJ Q 2 q t 

Thus the loading stiffness matrix for the last element becomes 


(3.31) 


K e L = 


<9f! 


dq t 

_ tCtT ^Qi 

L aqt 

= Nl Q 2 


jyjT dQl _ , tCtT 


3q t 


Qt + Nl Q 2 


(3.32) 
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Since the load only acts on the last finite element, the global dimensionless loading 
stiffness matrix is then defined as 


K l = 


0 0 
o K£ 


(3.33) 


3.2.5 Tangent Stiffness Matrix 

From Eq. (2.99), the internal virtual work is expressed as 

8W int = JJse^TdQ (3.34) 

n 

where £j are the generalized dimensionless strains given by Eq. (2.56) and Tj the 
generalized dimensionless stresses given by Eq. (2.55). Using the definition of virtual 
strains given by Eq. (3.13), Eq. (3.34) becomes 

SW M = StfJJ^Tdn (3.38) 


where f; e nt is the internal force vector. 


The dimensionless tangent stiffness matrix is given by 

-T\ r r Qe T dT 


K e = 


^_nt 

dq, 


:-II^) TdCl+ U^^ dCl (3 ' 36) 


n 


| 1 = B sd 

dq t 


and — — - D c — — = D c B Bd 

9q t dq t 


Note that 


(3.37) 
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Then the tangent stiffness matrix becomes 


K e 


- //^ Tdfi+ // B " D 4 


dCl 


(3.38) 


-//§ 


II 


D c edn+ / / BL D c Bsd dQ, 


Recalling that B s( j = Bl + Bnli where Bl is constant but Bnl depends linearly 
on qt , gives the well known decomposition of the element tangent stiffness matrix 


K e = K^j + K|, + Kg 


(3.39) 


where K^, Kf,, and denote the dimensionless element linear, initial-displacement, 
and geometric stiffness matrices, respectively. These are given by 

r\ rl/2 

Ky = / / Bl D c B l dy dx (3.40a) 

J0 J- 1/2 

r 1 z'i/2 

K e D = / / {b£d c B nl + B£ l D c B l (3.40b) 

Jo J- 1/2 L 

Dc dy dx 

k — (3 - 40c) 

where the width and the length of the beam element are nondimensionalized to 
unity, and D c is the dimensionless equivalent bending-stiffness matrix. Also, note 
that <9B^ L /<9q t is independent of displacements. 

Let us define the stiffness contribution of the material stiffness as the elastic 
stiffness matrix (Doyle, 2001): 


Ke = Km + Kd 


(3.41) 
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The tangent stiffness matrix is symmetric for conservative loading. However, 
for nonconservative loading it is unsymmetrical, the reason being that the loading 
stiffness is added. The above formulation has also been derived using indicial 
notation and it is presented in Appendix C. The convenience of indicial notation 
is that it helps the programming. 

3.2.6 Interior Node Condensation 

The interior nodes do not connect with the adjoining elements in the assemblage. 
Therefore, such internal degrees of freedom are condensed at the tangent stiffness 
level and thus the stiffness matrix, mass matrix, and nodal load vector are only 
expressed in terms of the corresponding exterior node displacements of the element. 

In order to proceed, the displacement vector is rearranged by partitioning the 
relevant terms corresponding to external and internal degrees of freedom as follows: 

qt = {qi.q2} T ( 3 - 42 ) 

where q x are the exterior nodes to be kept and q 2 are the interior nodes to be 
condensed: 

qj = {ui,Vi,Wi,Ti,<t>i,(3i,U2,V2,W2,T2,<p2,P2} T (3.43) 

q 2 = {u 3 ,u 3 ,i>4,u;3,^4,'r3,r4,^3,/%} T ( 3 - 44 ) 

Now the load vector can be partitioned in a similar way: 

F={F 1 ,F 2 } T (3-45) 

The element stiffness matrix and mass matrix are also rearranged and partitioned 


3.2. DISCRETIZATION OF THE GENERALIZED PVW 


81 


according to Eq. (3.42) as follows: 


K 

M 


Ku 

k 12 

k 21 

k 22 

Mu 

Mi 2 

M 2 i 

M 22 


(3.46) 

(3.47) 


The Irons-Guyan Reduction (Irons, 1963; Irons, 1965; Guyan, 1965) for a system 
of equations is described in Appendix D. Thus the condensed stiffness matrix, 
condensed mass matrix, and condensed load vector are obtained as follows: 


k r 

= Ku — K 12 Kj 2 K 21 

(3.48) 

m r 

= Mu — Ki 2 Kj 2 M! 2 i — N4i 2 K2 2 K 2 i 

(3.49) 


+ Ki 2 K 2 2 M 22 K 2 2 K 2 i 


f r 

= Fi — K 12 K 22 F 2 

(3.50) 


The interior node displacements can be obtained using 

cfc = -K 12 K 22 2 q, (3.51) 

3.2.7 Equations of Motion 

Now we assemble the elemental matrices and vectors, and substitute the global 
matrices and vectors into Eq. (3.15). The GPVW is now a function of the gener- 
alized displacements q t , which are independent of each other. In other words, the 
virtual displacements <5q t are entirely arbitrary. Thus it follows that the integral 
can only be zero for all <5qt if and only if the coefficients of <5qt are identically zero 
(Meirovitch, 1997). Thus the equations of motion for the global system are given 
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by 


Eyyk 2 h 2 


2 1 

h 2 “ UH ' E„bh 2 E m i? M Z 


f- ft ft _ Pt OJ- p 

e F = ^K E q t + -§K G q t --^^K L q t + — ^ 


M — T q t (3.52) 


J yy ‘ 


where F is the global external force vector, K E the global elastic stiffness matrix, 
Kq the global geometric stiffness matrix, Kj, the global loading stiffness matrix, 
and M the global mass matrix. Assuming the length of each element is the same, 
the above equation can be expressed in terms of the total length of the beam as 
follows: 

(3.53) 


n e i 


where n e i is the number of elements used. Also, recall that r\ = £ e /h. If we plug 
the above expression into Eq. (3.52) and multiply the entire equation by we get 


= n h r i Ke qt + r i K g qt — Kl q t + M ^ ^ ( 3 - 54 ) 


E, 


yy 


Eyy h Tl el U) 


Let us define the nondimensional eigenfrequencies and buckling loads as 


ui‘ 


u = 


p? 


Eyy h 2 


(3.55) 


and the axial load as 


P = 


Pn 2 el r\ PL 2 


Eyy b Eyy b h 


(3.56) 


where £ is the total length of the beam. Let us redefine the dimensionless vectors 
and matrices as follows: 


P = ^Il f K e = r\ K e m4m 
Eyy Kl 


(3.57) 


3.3. VIBRATION AND STABILITY ANALYSIS 


83 


For the case of an axially compressed beam, the pre-equilibrium conditions yield 
the following initial stress (per unit length): 

T = {-N xx , 0,0,0, 0} T (3.58) 


Thus the geometric matrix can be written in terms of the buckling load as follows: 

(3.59) 

Liyy & 

Thus the equation of motion becomes 


77 Pn 2 el rj 

Kg = --^t Kg 


- 2 f •• 


F = K e q t — PKg q t -FK L q t +^M — q t 


(3.60) 


u 


3.3 Vibration and Stability Analysis 


In general, for the analysis of aircraft structures, the study regarding the stabil- 
ity of these structures cannot be ignored. The reason is that structures such as 
wings undergo large deflections and are, in general, subject to conservative as well 
as nonconservative forces. In this section, we present the equations for the free 
vibration and stability analysis of laminated structures subject to conservative and 
nonconservative loading using the finite element method. 

3.3.1 Equation of Motion About the Equilibrium State 


A structure is stable at an equilibrium state if for every small disturbance of the 
system the response remains small. Thus when studying the stability and vibration 
about the equilibrium state, the generalized displacements are perturbed by an 
infinitesimal displacement: 


qt = qo + qi 


(3.61) 
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where qo represents the displacements in the equilibrium state, and qi the infinites- 
imal perturbation. Eq. (3.61) is substituted into the equation of motion, Eq. (3.60). 
Note that the mass matrix, linear stiffness matrix, geometric stiffness matrix, and 
loading stiffness matrix are all independent of displacements. However, the initial 
displacement matrix is the only matrix dependent on the displacements. Moreover, 
this matrix depends on the linear and nonlinear strain-displacement matrices. 

Since the nonlinear strain-displacement matrix is linearly dependent on dis- 
placements, the nonlinear strain-displacement matrix can be expressed as 

Bnl — B No + Bnx (3.62) 

On the other hand, the linear strain-displacement matrix is independent of displace- 
ments. As a consequence, the initial displacement stiffness matrix is expressed as 

Kd = JJ |Bl D c [B No -I- B Nl ] -I- [Bn 0 + B5J D c B l 

n 

+ [Bj[j 0 + Bjl^] D c [Bn 0 + Bnj] 

and expanding the above we get 

K d = JJ {B£ D c B No + Bj, 0 D c B L + B£ 0 D c B No } dCl 

« , 

s 11 v 

Kd 0 

+ JJ |Bl Dc Bnj + B^j Dc Bl + B^ D c Bn 0 4- Bn 0 D c Bnj j dCl 
n 

v ,i — v 1 

K Dl 

+ JJ {b^ D c B Ni | dCl 

Cl 
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Now the equation of motion can be written as 

Km qt + Kd qt — P Kg qt — P Kl qt + ^ M — qt — F — 0 

U) 

Km (qo + qi) - -PKg (qo + qi) ~ P k l (qo + qi) 

+ [K Do + K Dl + K Da ] (qo + qi) + w 2 M (q 0 + qi) ^ ~ F = 0 

and by expanding we get 

K m qo + K Do qo - P Kg qo - P Kl qo + M — j qo - F 

Lu 

v — — — ; 

vanishes because of equilibrium 

+ K M qi + K Do qi + K Dl q 0 - P K G qi - P K L qi + w 2 M -j qi 
+ Kd! qi + Kd 2 qo + Kd 3 qi = 0 

Assuming that a trivial, rotationless, equilibrium state exists, the out-of-plane dis- 
placements are zero. Note that the initial displacement matrix is only a function 
of the out-of-plane displacements. Under this assumption the equation of motion 
reduces to 

K m qi — P K g qi — P K L qi + 0? M \ qi + K Dl qi + K Da q x - 0 (3.63) 

UJ S v 

s v* x Nonlinear 

Linear 

where K G is the initial stress matrix and K L is the loading matrix. The elastic 
stiffness matrix becomes 


Kg = Km + Kdi H" Kd 3 (3.64) 

Moreover, for a wide range of problems it is of great interest to study the linearized 
buckling and vibration problem. For the linearized case we neglect the nonlinear 
terms in Eq. (3.63). In other words, we assume that the change in geometry prior 
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to the first critical state can be neglected and the initial displacement matrices are 
neglected, i.e., 

K e = K m (3-65) 

Further assuming a harmonic motion with a frequency of u, 

qi = <!>e iu \ (3.66) 

the equations of motion to study the linearized stability and vibration about an 
equilibrium state can be expressed as 

K m <f> — P K g <}>- P Kl <t> - w 2 M <j> = 0 (3.67) 

3.3.2 Free Vibration Analysis 

The free vibrations about the trivial equilibrium state are studied by taking P = 0. 
Moreover, recall that we assumed a harmonic response, thus equations of motion 
result in an eigenvalue problem: 

[K m - C? M] <t> = 0 (3.68) 

where w 2 is the dimensionless eigenfrequency, and 4> the corresponding dimension- 
less right eigenvector. 

3.3.3 Linearized Buckling Analysis 

In general, for all types of forces the dynamic analysis will always predict the buck- 
ling load. The dynamic criterion considers small oscillations about the equilibrium 
position, and reduces the stability problem to that of solving an eigenvalue problem 
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to determine natural frequencies (ui) and associated eigenmodes, i.e., 


K m - FK g - FK l - cD 2 M 


0 = 0 


(3.69) 


Note that the frequencies are generally complex numbers. The system’s dynamic 
stability depends on the values of u> 2 : 


(a) u> 2 > 0, and purely real: The system is dynamically stable 

(b) u> 2 = 0: The system is dynamically critical 

(c) Unstable otherwise. 

The dynamic criterion can be applied to both conservative and nonconservative 
systems. 


Conservative Loading 

For the case of conservative loading, a load potential function usually exists and 
K l = 0. Thus we can use the static criterion, which looks at admissible static 
perturbations of an equilibrium state. This results in the following linearized eigen- 
value problem: 

Km — P Kg j 0 = 0 (3.70) 

Here we find a P such that the lowest eigenvalue of the tangent stiffness matrix is 
zero. In light of the dynamic criterion, the buckling load is found such the system 
is dynamically critical. Thus the eigenvalue problem to be solve is 

K m -PK g -w 2 m]0 = O (3.71) 

where we find a P such that u) 2 = 0. 
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Figure 3.6: Subtangential Loading. 


Subtangential Loading 

For the generally nonconservative loading, consider a cantilevered laminated beam 
subject to a combined conservative load F con and a tangential follower force F nc at 
the tip, as shown in Fig. 3.6. The former may be thought of as a dead load, while 
the latter can be considered as a typical follower force. When working with sub- 
tangential loading, where the conservative load is combined with a nonconservative 
one, it is convenient to define a nonconservativeness parameter (Langthjem, 2000): 

r) - —r~ Ft Fon + Ac (3.72) 

Ft 

where 77 G R and is the nonconservativeness loading parameter, P nc the dimension- 
less nonconservative force, P co n the dimensionless conservative force, and Ft the 
magnitude of the sum of the conservative and nonconservative loads. Note that the 
value 77 = 0 corresponds to a pure dead load (conservative load), rj / 1 corresponds 
to a generally nonconservative load, and 77 = 1 corresponds to a purely tangential 
follower load (nonconservative load). 

In general, it is sometimes convenient to define a loading matrix, L, as the 
contribution of the geometric stiffness and the loading stiffness correction: 


L = Kg - v K l 


(3.73) 
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The loading matrix does not depend on the material nor on lamina characteristics 
such as ply angle and ply thickness. It is only a function of length, width, and the 
nonconservativeness parameter. Moreover, it should be noted that some authors 
may define the loading matrix as L = Kg + t/Kl- The reason for the sign difference 
is that they define the loading stiffness, given by Eq. (3.23), as the negative of 
the partial derivatives of the nonconservative loading vector with respect to the 
generalized displacements. Although both expressions are equivalent, in this work 
we will use the notation as given by Eq. (3.73). 

Thus for the case of subtangential loading, the linearized eigenvalue problem 
can be expressed as 

Km-PL-^mJ^O (3.74) 

Here the buckling load can be by divergence and/or flutter. We compute all the 
eigenvalues of Eq. (3.74) at each load step in correspondence to the deformed 
configuration. Note that the total tangent stiffness matrix at each load step may 
be unsymmetric for every t] > 0. Thus the stability conditions can be written as 

(a) if all Qj 2 n are real and positive numbers, then the equilibrium state is stable. 
The critical load is the value of the load for which the smallest eigenfrequency 
becomes zero, iLf = 0; 

(b) if all t2>n are real numbers, and CJ\ = 0 becomes zero, then the equilibrium state 
becomes unstable and the stability transition occurs via divergence instability. 
The critical load is the value of the load for which the smallest eigenfrequency 
becomes zero, Cj\ = 0; 

(c) if at least one pair of the eigenvalues wfn-i A becomes complex conjugate, 
then the equilibrium state is unstable and the stability transition occurs via 
flutter instability. The critical load is the value of the load for which two 
eigenfrequencies approach each other until they coalesce, wln-i = &\ n . 
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Thus, we find P such that Qj\ n _ 1 = u> 2 n or Cb\ — 0. 


3.4 Summary 


In this chapter, we have discretized the generalized principle of virtual work. By 
doing do, we developed the element and global finite element matrices and vec- 
tors. These were obtained through exact integration using MATHEMATICA as 
the symbolic manipulator. 

We also derived the equations for both free vibrations and the stability of the 
equilibrium state for laminated beams subject to subtangentially loaded beams. 


Chapter 4 
Deterministic Finite Element 
Analysis of Laminated Beams 


In general, for the analysis of aircraft structures, nonlinear analysis should be 
considered because structural components such as wings undergo large deflections. 
The equations of motion for the study of large deflections were derived in Chapter 
3 using the weak form of the generalized principle of virtual work. 

For a conservative system, the only possible initial instability is of divergence 
type. For a nonconservative system, however, instability can be by divergence, 
flutter, or both, depending on the amount of nonconservativeness. In this chapter, 
we present various results for the static and free vibration analysis of various cases 
present in the literature. Then we use the dynamic method to study the stability 
of conservative and nonconservative structural systems. 


4.1 Case properties and definitions 

The purpose of this section is to explain how the results are obtained. First we dis- 
cuss how the results are nondimensionalized, then the various boundary conditions 
used. 


91 


4.1. CASE PROPERTIES AND DEFINITIONS 


92 


4.1.1 Nondimensionalization 


In the literature, researchers nondimensionalize eigenvalues and buckling loads in 
two different forms. Thus throughout this dissertation the natural frequencies are 
obtained in their nondimensionalized form in either of the following forms 

Eyy b 0 kg 

and the critical loads as 

e 

p — p 

n n p h h 3 

f-Jyy "o ' *'o 

Note that the and P n are equivalent to those given in the previous chapter. 

Let the dimensionless left and right eigenvectors of the fc th mode be defined by 
ipy and <^>k) respectively. In order to normalize the left and right eigenvectors for 
nonconservative systems, we need two independent criteria. Thus let us normalize 
the eigenvectors such that 



12 Jo t 


Eyyb'hl 


(4.1) 


Pn = Pn 


12 r 


Eyy b Q h Q 


(4.2) 


{0k} [M] {0k} — 1 and {0k} n tt> nonzero element {0k} rl <-li nonzero ele:nent 

for a selected value n. In other words, we know that we can multiply the eigenvectors 
by any arbitrary scalar: 


{0k} = a {0k} {0k} = b {0k} 


Then the second condition suggests 
a0 k (n) = 60k(n) 


_ . 0k (n) 

a 0k («) 


(4.3) 
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and the first one suggests 


ab {0k} T |M] {0k} = 1 


{0 k } T |M]{0 k } 


Using Eqs. (4.3) and (4.4), we get the normalization constants a and b : 


(4.4) 


a = 


l<h kW 1 

^k(n) {t/> k } T [M] {<ftj 


b = 


/ V>k(n) 1 

<t> k(«) {V> k } T [M] {</) k } 


(4.5) 


For conservative problems, the left and right eigenvectors are equal. Thus only one 
criterion is needed and in fact 


a = b = 


I 1 
{<f> k } T [M] {0 k } 


(4.6) 


4.1.2 Boundary Conditions 

For the types of problems treated in this dissertation, we consider three sets of 
boundary conditions. These sets of boundary conditions are shown in Fig. 4.1 and 
are given as: hinged-hinged, clamped-free, and clamped-clamped. In all cases the 
load is applied at the tip of the beam {x — E). The consequences of these boundary 
conditions on the generalized midplane displacements are: 

Clamped-Free: u = v = w = r = (j) = /3 = 0(x = 0) 


u = v = w = 
v = w = 


t = <p — P = 0 (x = 0) 

T — (f> — P — 0 (X = i) 


Clamped-Clamped: 
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Figure 4.1: Boundary 


DEFINITIONS H 


Hinged-Hinged 



Clamped-Free 


Clamped-Clamped 



used for the analysis of laminated beams. 
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Hinged-Hinged: 


u = v = w = t = 0 (x = 0 ) 
v = w = t = 0 (x = £) 


For nonconservative problems, only cantilever (clamped-free) laminated beams were 
considered. 


4.1.3 Computer Program: NLbeam21.f 

In this work, the symbolic program MATHEMATICA is an interface to FORTRAN 
90. The MATHEMATICA file NLbeamSl.nb is used to calculate all the element 
matrices and vectors symbolically and printed in Fortran-form. Once these matrices 
are available they are then inserted into the FORTRAN program NLbeam21.f. 

The FORTRAN program NLbeamSl.f is written entirely in double precision 
and works as the pre-processor, processor, and post-processor. The code reads the 
input file NLbeam21.dat, which provides all the necessary information to perform 
the wanted analysis. The program NLbeam21.f uses LAPACK 90 libraries 1 to 
solve the eigenvalue problem. A copy of the code can be obtained by contacting 
the author or the Aerospace and Ocean Engineering Department here at Virginia 
Tech. 


4.1.4 Stability Analysis Using ABAQUS 

The buckling and flutter loads obtained with our computer program are verified 
with those using the finite element package ABAQUS. The preprocessing is done in 


1 These libraries are available online: http://www.fortran.com 
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PATRAN and the output is fed into ABAQUS for processing and post- processing. 

We had no problem in finding the buckling load in ABAQUS. On the other 
hand, the flutter load was challenging because ABAQUS (i) cannot solve for flutter 
load, (ii) its eigensolver has problems when frequencies are very close to each other 
(ABAQUS User’s Manual, 1998), and (iii) requires a lot of elements to achieve 
convergence. 

In order to approximate the critical load in ABAQUS, we found an alternative 
way of doing so. First, we apply a small eccentricity to the applied axial load by 
adding an offset to the midsurface. This can be achieved by adding the OFFSET 
command in the ABAQUS’ input file as follows: 

♦SHELL SECTION, COMPOSITE, ELSET=LAMBEAM , 0FFSET=-0.1 


In order to model subtangential forces, say r] = 0.6, we add the following lines to 
the input file: 


♦NSET, NSET=FUERZA , GENERATE 
151, 156, 1 

♦STEP, UNSYMM=YES , NLGEOM , INC=10000 
♦STATIC 

1 . 0 , 100 , , 1.0 

♦CL0AD, FOLLOWER 

FUERZA, 1, -60.0 

♦CL0AD 

FUERZA, 1, -40.0 

♦RESTART, WRITE, FREQ=2 
♦NODE FILE 
U, CF 

♦ENERGY FILE 
♦ENERGY PRINT 
♦END STEP 
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Wtip 

0.003 

0.002 
0.001 
0.000 
- 0.001 

- 0.002 

- 0.003 
- 0.004 
- 0.005 
- 0.006 
- 0.007 

Figure 4.2: Approximating t.lit* flutter load in ABAQIJS. 

Finally, we perform a “post-buckling” type analysis to determine the flutter load. 
Figure 4.2 illustrates this method. We know that at the onset of Hut ter, the trans- 
verse displacement at the tip of the cantilevered beam is minimum. After we have 
reached flutter and we increase the load, the beam will shift modes. In other words, 
the tip displacement will increase. This will be the flutter load. 

It should be noted that ABAQUS does not have a laminated beam element, 
and the only possibility is to use a shell element, such as S4R, to model the can- 
tilevered laminated beam. Some advantages of our program are that we can plot 
the fundamental characteristic curves and solve for the divergence or flutter loads 
without problems because the eigensolver built, in LAPAC'K has no problems when 
two frequencies are close to each other (even when they are identical), and only five* 
elements are required to achieve convergence. 
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4.2 Case I: Conservative 


To study the stability of conservative systems using the dynamic criterion, the 
following eigenvalue problem is solved: 



- P K g — clrM 


<t> = 0 


We find a P n such that Cj\ = 0. In the free vibration analysis we find a for 
P = 0. 


4.2.1 Isotropic Beams 


To validate the present shear deformable laminated beam element, we studied the 
cases for isotropic beams and compared the results to those by Euler (classical 
buckling loads). The laminated beam element was treated as an isotropic material 
with the following mechanical properties: 


V xy — ^xz V yz 0 


Only one ply 
properties: 


Et 


= 1.00 


G 


xy 


Gz 

E„ 


v z 


E„ 


0.50 


J yy ^vv vy 

was considered for the isotropic case with the 


following geometric 


— = 0.5 b a = 0.1 ft 


For the isotropic case, any value of the ply angle produces the same result. All 
results are compared to those by an analytical solution (Meirovitch, 1997; Jones, 
2001 ). 
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Cantilevered Beams 

The first two dimensionless natural frequencies for the isotropic cantilevered beam 
are 


Analytical: u)i = 3.516 u> 2 — 22.03 

Present: u)\ = 3.516 £>2 — 22.03 

The first two dimensionless critical loads for the isotropic cantilevered beam are 

Analytical: Pi = 2.467 P 2 = 22.21 

Present: Pi — 2.467 P 2 — 22.22 

Figure 4.3 shows that the mode shapes are in perfect agreement with those given 
in the literature. 

Clamped-Clamped Beams 

The first two dimensionless natural frequencies for the isotropic clamped-clamped 
beam are 


Analytical: £b\ = 22.37 £2 — 61.70 

Present: u) 1 = 22.35 ^2 — 61.70 

The first two dimensionless critical loads for the isotropic clamped-clamped beam 
are 

Analytical: Pi = 39.48 P 2 = 80.73 

Present: Pi — 39.52 P 2 = 81.32 
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Figure 4.4 shows that the mode shapes are in very good agreement with those given 
in the literature. 

Simply-Supported Beam 

The first two dimensionless natural frequencies for the isotropic simply-supported 
beam are 


Analytical: uJ\ — 9.870 <^2 — 39.48 

Present: = 9.870 = 39.52 

The first two dimensionless critical loads for the isotropic simply-supported beam 
are: 

Analytical: Pi = 9.870 P 2 = 39.48 

Present: Pi = 9.870 P 2 = 39.52 

Figure 4.5 shows that the mode shapes are in very good agreement with those given 
in the literature. 

Results for both dimensionless natural frequencies and buckling loads are in 
good agreement with those of the exact solutions for all three boundary conditions 
mentioned above. 

4.2.2 Laminated Beams 

Next, the results for various laminated beams were compared with those by Maiti 
and Sinha (1994), and Reddy (1997). The mechanical properties for laminated 
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■ 1st Mode: 5 elements 
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2nd Mode: 5 elements 


(a) Vibration mode shapes 
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(h) Buckling mode* shapes 


Figure 4.3: Mode shapes for an isotropic beam with clamped- free boundary conditions 
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(a) Vibration mod** shapes 



(b) Buckling mode shapes 

Figure 4.4: Mode shapes for ail isotropic beam with clamped-clainped boundary condi- 
tions. 



(a) Vibration modi* shapes 



(b) Buckling mode shapes 


Figure 4.5: Mode shapes for an isotropic beam with hinged-hinged boundary conditions 
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beams as given by Maiti and Sinha (1994) are 

— = 13.70880 ^ = 0.54710 ^ = 0.45679 

Eyy Eyy Eyy 

with 

kg 

Eyy = 9.42512 x 10 9 Pa p = 1550.0666 

ww m 

and the geometric properties are 

^ = 0.3175 b a - 0.01 m 

b 0 

and the slender ratio t/h Q and ply angles 0 l are individually specified for each 
problem. Here h Q represents the total thickness of the beam and each ply is assumed 
to have the same thickness. Thus the thickness of each ply is calculated as h n = 
h 0 /N\ am , where iV lam is the total number of plies. 


a 

E, 


v — = 0.26964 


yy 


v xy = 0.30 


The mechanical properties for laminated beams as given by Reddy (1997) are 


E X x 

Eyy 


= 25.00 


^ = 0.50 ^ = 0.50 

Eyy Eyy 



0.20 


with 


E y y = 1.9584 x 10 8 psf 


p = 0.250387 


slugs 

ft 3 


v xy = 0.25 


and the geometric properties are 

— = 0.50 b a = 0.1 ft 

b 0 

and the slender ratio t/h 0 and ply angles 0 { are individually specified for each 
problem. Thus the thickness of each ply is calculated as h n = h 0 /N \ am , where N] &m 
is the total number of plies. 
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The ply and the beam’s mechanical and geometrical properties are assumed to 
be uniform throughout the beam. The properties are given as a footnote to the 
table in which results for a given example are presented. All results are obtained 
with five elements. 

We proceed to calculate the dimensionless natural frequencies for various unidi- 
rectional laminated beams and all three boundary conditions previously mentioned. 
Table 4.1 contains the results with warping and without warping. Results shows 
that the dimensionless natural frequencies are in very good agreement with the re- 
sults obtained by Maiti and Sinha (1994). In fact, although they used a nine-noded 
rectangular isoparametric element, our results using a one-dimensional FSDT beam 
theory compare well with theirs. 

Also, results show the importance of including the weirping function in a lam- 
inated beam element. When warping is neglected, the bending term ( b Q D^) gets 
modified to (b a Dee + A 55 bl/12), leading to very high frequencies. The proposed 
warping function removes this addition to D 66 , and thus leads to very good results. 
The small discrepancy has to do with the fact that we are modeling the plate strip 
using a beam element, whereas Maiti and Sinha modeled it using a plate element. 

For all cases shown in Table 4.1, it can be seen that as we increase the ply 
angle from 0° to 90°, the dimensionless natural frequency decreases. Thus the fiber 
orientation plays an important role in the design of laminated structures. One must 
be aware that the fundamental natural frequency may decrease or increase with a 
change in the ply angle. 

In Tables 4.2, 4.3, and 4.4 we compare various cases of laminated beams studied 
by Reddy (1997). The dimensionless fundamental natural frequencies and critical 
buckling loads are in good agreement with those obtained by Reddy (1997). The 
results for the simply-supported laminated beam seem to compare not that well 
with those by Reddy (1997). One reason could be that Reddy (1997) does not take 
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into account the various coupling effects as is the case for the present study. On 
the other hand, our results are in good agreement with those presented by Maiti 
and Sinha (1994), see Tables 4.5 and 4.6. 

Table 4.5 shows results for various symmetrically and unsymmetrically lami- 
nated beams. Once again, results are in very good agreement with those by Maiti 
and Sinha (1994). The results show that as the length-to-thickness ratio increases, 
the dimensionless natural frequency approaches an asymptotical value. A com- 
parison of Tables 4.1 and 4.5 shows that when unidirectional laminated beams 
are sandwiched with two plies of 0°, while keeping the beam s material and ge- 
ometric properties the same, the natural frequency increases when compared to 
undirectional laminated beams. Moreover, the frequencies are very close to those 
by obtained with a 0° laminate. 

We also present the buckling loads for each of the cases presented by Maiti 
and Sinha (1994). These results are tabulated in Table 4.7. As in the case of free 
vibrations, the results show that the warping effects lead to lower dimensionless 
buckling loads. This is especially the case for unidirectional plies of 6 = 30° and 
45°. Also, the buckling load increases as the unidirectional plies are sandwiched 
with two plies of 0°, while keeping the the beam’s mechanical properties the same. 

Since results presented in Tables 4.1, 4.5, and 4.7 show that warping is extremely 
important for a unidirectional beam with a ply of 30° layout, we also checked the 
first two mode shapes related to this case. Figure 4.6 shows that the fundamental 
vibration and buckling modes are not affected by warping, although the natural 
frequencies and buckling loads are. However, the second buckling mode is slightly 
affected. This shows that warping can affect higher buckling loads as well as the 
mode shapes. 

All these results were obtained using our dimensionless finite element formu- 
lation. In the dimensional model, the stiffness matrix may have coefficients with 
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(a) Vibration mode shapes 



xA 


(b) Buckling mode shapes 

Figure 4.6: Mode shapes for a cantilevered unidirectional laminated beam (0 = 30°). 
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large differences in the order of magnitude from than others within the matrix. The 
same may be true for the mass matrix. An advantage of the dimensionless model 
over the dimensional one is that the elements of the stiffness matrix are of the same 
order of magnitude. The same is also true for the mass matrix. The results for 
the cases presented by Goyal and Kapania (2002) using a dimensional finite ele- 
ment formulation were in perfect agreement with those presented here when using 
a dimensionless finite element analysis. 

We also studied the case when a unidirectional laminated beam is subject to 
a pure dead load weight (77 = 0) with warping included (a = 1). The mechanical 
properties for the single-ply laminated beam, as given by Xiong and Wang (1987), 
are 

— = 3.4339623 ^ = 0.4462264 ^ = 0.4462264 ^ = 0.3707547 

Eyy Eyy Eyy Eyy 

with 

Eyy = 0.10388 x 10 9 Pa p= 1860.0 u xy = 0.33 

w m 

and the geometric properties are 

— = 0.25 1- = 10.0 b Q = 0.02 m 

b a h 0 

and the ply angle 6 takes values of 0°, 15°, 30°, 60°, and 90°. Here h Q represents 
the total thickness of the beam. 

Figure 4.7 shows the characteristic curves for the fundamental buckling load 
for unidirectional cantilevered laminated beams subject to a conservative load- 
ing. Results show very good agreement with those obtained using ABAQUS. Plots 
show that the ply orientations play a significant role in the stability of unidirec- 
tional laminated beams. For the conservative case, the fundamental dimensionless 
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Figure 4.7: Characteristic curves for tin 1 fundamental buckling loud P for unidirectional 
cantilevered laminated beams subject to a conservative loading. 


natural frequency Cj\ and the first buckling load I\ are bounded by those of the 
unidirectional laminated beam with orientations of () u and 90°. 

Since it w<us shown that warping is most important for unidirectional laminated 
beams with a play angle of 30°, we proceeded to compare the buckling mode shape 
obtained with our laminated beam (dement and with that using ABAQUS. Figure 
4.8 shows that when using only five laminated beam (dements are sufficient to 
capture the same effect obtained by ABAQUS. Results in Fig. 4.8(b) where obtained 
using MATHEMATICS. 
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Figure 4.8: Twisting effect, in the buckling mode captured using ABAQUS and present 
formulation. 
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4.3 Case II: Generally Nonconservative 

Since it was shown in the previous section that warping effects are important in the 
analysis of laminated composite beams, hereon the analysis will include warping. 
In other words, we take the warping constant to have the following value: 

a = 1 


When the beam is subject to subtangential loading, the linearized eigenvalue 
problem can be expressed as 

K m -PL-u 2 m]^ = 0 


where 

L = K g - V K l 

and 7 } is the nonconservativeness loading parameter. Here we find P„ such that 
(b\ n _i = Cj\ n or C? n = 0. In the free vibration analysis we find for P = 0. 

4.3.1 Isotropic Beams 

The geometric and material properties for isotropic beams are the same as those 
used for the conservative case. Here we consider the case of a generally noncon- 
servative system. The load becomes nonconservative as the nonconservativeness 
loading parameter r? takes values other than zero. We plotted several character- 
istic curves corresponding to different values of r). Results are shown in Fig. 4.9 
and are in perfect agreement with those presented by Gasparini et al. (1995). Al- 
though Gasparini et al. (1995) used twenty finite elements for their analysis, here 
five elements were sufficient for convergence. 
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Figure 4.9 is composed of two figures. Figure 4.9(a) shows the complete plot 
of the characteristic curves, whereas Fig. 4.9(b) shows the regions of interest in an 
expanded form. As seen in these figures, the buckling load decreases as the loading 
parameter decreases. As Fig. 4.9 shows, the fundamental characteristic curves for 
the various values of 77 intersect at the flutter load when 77 = 0.5. The dimensionless 
flutter load at rj = 0.5 is P = 16.2. 

The value of the nonconservative parameter 77 for a particular problem is very 
important in design. The reason is that by knowing how nonconservative the load 
is, one can predict if the instability will occur due to divergence or flutter. In 
the case of cantilevered isotropic beams with constant cross sectional material and 
geometric properties, and in the absence of damping, all values of rj < 0.5 will lead 
to a divergence type instability; and to a flutter type instability otherwise. 

The flutter load and fundamental frequency modes for the case of purely non- 
conservative load are shown in Fig. 4.10. 


4.3.2 Laminated Beams 

In the literature, little attention has been paid to the effect of the nonconservative- 
ness loading parameter 77 on the stability of laminated beams. Xiong and Wang 
(1987) studied the stability of a Beck-type laminated column using an analytical 
model, which is a column subject to a purely tangential load. Thus here we first 
validate our finite element model by comparing our results to the analytical solution 
of Xiong and Wang (1987). 

The mechanical properties for the single-ply laminated beam as given by Xiong 
and Wang (1987) are 

— = 3.4339623 — = 0.4462264 ^ = 0.4462264 ^ = 0.3707547 

Eyy Eyy Eyy Eyy 
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with 

kfr 

E vv = 0.10388 x 10 9 Pa p= 1860.0 — v xy = 0.33 

vv m 

and the geometric properties are 

— = 0.25 4- = 10.0 b a = 0.02 m 

b Q h 0 

and the ply angle 6 takes values of 0°, 15°, 30°, 60°, and 90°. Here K represents 
the total thickness of the beam. 

We considered the case when the beam is subject to a pure tangential follower 
load (7] = 1). Figure 4.11 shows the fundamental characteristic curves for unidirec- 
tional cantilevered laminated beams subjected to a tangential follower load. As in 
the case of conservative loading, discussed in the previous section, plots show that 
the ply orientations play a significant role in the stability of unidirectional laminated 
beams. Results for the nonconservative case are in good agreement when compared 
to those by Xiong and Wang (1987) and ABAQUS. Furthermore, the buckling and 
vibration modes are shown in Fig. 4.14. As results show, the buckling and vibration 
modes remain unchanged as we vary the ply angle. Thus, although the critical load 
at flutter may vary with the ply orientation the mode shape remains unchanged. 

Now we proceed to study the effect of subtangential loading on unidirectional 
laminated beams (0 < rj < 1). Figure 4.12 shows results for a nonconservativeness 
loading parameter of 0.5 and Fig. 4.13 shows results for a nonconservativeness 
loading parameter of 0.35. When r] = 0.5, the instability occurs by divergence for 
any unidirectional laminated beam, as was the case for isotropic beams. It can 
also be observed that the critical load due to flutter is lower when compared to 
the purely nonconservative case (77 = 1). As was the case of isotropic beams, for 
values of 77 < 0.5 instability is governed by divergence in the case of unidirectional 

laminated beams. 
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Dimensionless P 


(a) C-liararlrristif nirvrn 



Dimensionless P 


(b) Region of intert;st 


Figure 4.11: Fundamental characteristic curves for unidirectional cantilevered laminated 
beams subject to a tangential follower load. The nondimensional load is P. 
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(h) Region of interest 

Figure 4.12: Fundamental characteristic curves for unidirectional cantilevered laminated 
beams subject to a subtangential follower load (tj - 0.5). The noudimensional load is P. 
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Figure 4.13: Fundamental characteristic curves for unidirectional cantilevered laminated 
beams subject to a subtangential follower load (// = 0.35). The nondiniensional load is 

P. 
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Figure 4.14: Modes at the onset of flutter for unidirectional cantilevered laminated beams 
subject to a purely tangential follower load. 
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Since it was shown that warping is most important for unidirectional laminated 
beams with a play angle of 30°, we proceeded to compare the flutter mode shape 
obtained with our laminated beam element and with that obtained using ABAQUS. 
Figure 4.15 shows that only five laminated beam element are sufficient to capture 
the same effect obtained by ABAQUS’ shell element S4R (with 125 finite elements). 
Results in Fig. 4.15(b) where obtained using MATHEMATIC A. 


4.4 Stability of Unsymmetrical Laminated Beams 


We also studied the stability of several cases of unsymmetrical laminated beams. 
The cases studied have a ply stacking sequence of [0°/#] where 6 takes values of 0 , 
15°, 30°, 60°, and 90°. We considered four cases: (i) purely dead load (*7 = 0), (ii) 
purely tangential follower load (ry = 1), (iii) subtangential load with 77 — 0.5, (iv) 
subtangential load with 77 = 0.35. For this analysis we used the properties as given 
by Xiong and Wang (1987) and compared our results with those obtained using 

ABAQUS. 

We first studied unsymmetrical laminated beams with a purely conservative 
loading. The results are shown in Fig. 4.16. Then we studied the case when the 
loading is generally nonconservative. The results are shown in Figs. 4.17, 4.18, and 
4.19. All results compare well with the ones obtained using ABAQUS. Results show 
similar behavior when compared to those of the unidirectional case. 
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Figure 4.15: Twisting effect, in the flutter mode captured using ABAQUS and present 
formulation. 
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Dimensionless P 

Figure 4.16: Characteristic* curves for the fundamental buckling load P tor unsvmmetrical 
cantilevered laminated beams subject to a conservative loading. 

4.5 Summary 

In this chapter, we have studied the deterministic stability and vibrational response 
of isotropic and laminated beams. The present twenty-one degrees ot freedom lam- 
inated beam element has been validated with those available in the literature. For 
the case of stability analysis, using the dynamic criterion, a very good agreement 
exists with the analytical and numerical analyses available in the literature. Re- 
sults have also shown the importance of ply orientation and warping effects when 
studying the stability of laminated beams. 

The analysis performed in this chapter is only valid for deterministic systems. 
However, in the presence of uncertainties a more rigorous method must be used. In 
the next chapter, we discuss the probability approach used in this dissertation. 
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Figure 4.17: Fundamental characteristic curves for unsymmetrical cantilevered laminated 
beams subject to a tangential follower load. The noudimensional load is P. 
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(a) Characteristic curves 



(h) Region of interest, 


60. Off 


60.00 


Figure 4.18: Fundamental characteristic curves for unsymmetrical cantilevered laminated 
beams subject to a subtangential follower load ( 7 / = 0.5). The nondimensional load is P. 
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Dimensionless P 
(h) Characteristic curves 



Dimensionless P 

(h) Region of interest 

Figure 4 . 19 : Fundamental characteristic curves for unsymmetrical cantilevered laminated 
beams subject to a subtangential follower load (?/ = 0.35). The nondimensional load is 

P. 
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Table 4.1: Dimensionless fundamental frequencies (u>) for unidirectional laminated beams 
with various boundary conditions (ct = 0, implies warping is not included, oc 1, implies 
war ping is included) 

0 tjho = 20 £/h 0 = 60 

Maiti & Present Maiti & Present 

Sinha (1994) a = 0 a = 1 Sinha (1994) a = 0 a = l 


Clamped-free 



0° II 

12.6450 

12.7908 

12.7908 

12.9940 

12.9925 

12.9940 


o 

o 

CO 

5.3374 

7.0673 

5.2928 

5.3926 

7.1125 

5.3113 


45° 

4.1627 

4.6571 

4.1489 

4.1815 

4.6716 

4.1592 

— 

05 

o 

0 

3.6981 

3.7705 

3.6930 

3.7061 

3.7793 

3.7020 


CO 

o 

o 

3.5103 

3.5068 

3.5068 

3.5154 

3.5150 

3.5150 

— 

Clamped- Clamped 







0° 

70.9066 

70.6606 

70.6606 

81.5849 

81.2112 

81.2112 



30° 

34.0387 

42.6153 

33.2170 

36.1359 

44.9922 

33.8174 


45° 

26.5628 

28.8494 

25.9482 

27.6692 

29.6489 

26.4301 


o 

O 

CO 

23.3448 

23.5125 

23.0653 

23.9469 

24.0046 

23.5178 

— 

90° 1 

21.9579 

21.8625 

21.8625 

22.4265 

22.3255 

22.3255 


Simply-Supported 







0° 1 

35.2637 

35.2325 

35.2325 

36.4265 

36.3930 

36.3930 


CO 

o 

o 

16.4237 

21.3061 

19.8201 

18.4462 

21.5981 

19.9389 


45° 

12.1022 

13.4558 

13.1341 

12.7466 

13.5444 

13.2020 


60° 

10.4288 

10.6078 

10.5673 

10.5659 

10.6602 

10.6176 


90° 

9.8248 

9.8159 

9.8159 

9.8735 

9.8645 

9.8645 

— 

E xx /L 

5™, = 13.7088, G xy /E yy = 

0.5471, G yz 

/Eyy = 0.269641 

> Gxz/Eyy 

= 0.45679, 


v xy = 0.30, h 0 /b 0 = 0.3175, E yy = 9.42512 x 10 9 Pa, b a = 0.01 m, p = 1550.0666 kg/m' 
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Table 4.2: Dimensionless fundamental frequencies (cl>) 
(A) for symmetrically laminated beams with t/h 0 — 
included; a = 1, implies warping is included) 


and dimensionless buckling loads 
10 (a = 0, implies warping is not 


e 

(deg) 


Reddy 

(1997) 


A 

Present Reddy 

(q = 0) (a = 1) I (1997) 


u 

Present 


(a = 0) (a = 1) 


Clamped- Free 


0 

4.576 

4.576 

90 

0.203 

0.203 

[0/90], 

3.922 

3.922 

[45/ - 45], 

0.355 

0.360 


Clamped- Clamped 


0 

27.656 

27.686 

90 

2.747 

2.755 

[0/90], 

20.800 

20.819 

[45/ - 45], 

4.802 

4.870 


Simply- Supported 


0 

13.768 

13.770 

90 

0.784 

0.784 

[0/90], 

11.179 

11.181 

[45/ - 45], 

1.369 

1.441 


4.576 

4.528 

4.560 

4.560 

0.203 

1.004 

1.002 

1.002 

3.922 

4.132 

4.178 

4.178 

0.355 

1.326 

1.333 

1.332 

27.686 

17.212 

17.215 

17.215 

2.755 

5.761 

5.764 

5.764 

20.819 

14.837 

14.839 

14.839 

4.815 

7.616 

7.666 

7.623 

13.770 

11.635 

11.636 

11.636 

0.784 

2.771 

2.771 

2.771 

11.181 

10.488 

10.488 

10.488 

1.437 

3.663 

3.758 

3.752 


E xx / E ] 


yy 


= 25, Gxy/Eyy = 0.5, G y Z / Eyy = 0.2, G XZ /Ey y = 0.5, V X y = 0.25, h 0 /b 0 = 0.5, 
Eyy = 1.9584 x 10 8 psf, b a = 0.1 ft, p = 0.250387 slugs / ft 
(j Eyy = 9.37687 x 10 9 Pa, b a = 0.03048 m,p= 129.0290 kg/m 3 ) 
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Table 4.3: Dimensionless fundamental frequencies (<i) and dimensionless buckling loads 
(A) for symmetrically laminated beams with i/h a = 20 (a = 0, implies warping is not 
includ ed; a = 1, implies warping is included) 

A £» 

6 Reddy Present Reddy Present 

(deg) 1 (1997) (a = 0} (a = 1) | (1997) (a = 0) (a = 1) 


Clamped-Free 


0 

4.987 

4.987 

4.987 

4.930 

4.931 

4.931 

90 

0.205 

0.205 

0.205 

1.012 

1.012 

1.012 

[0/90] s 

4.362 

4.362 

4.362 

4.594 

4.597 

4.597 

[45/ - 45], 

0.358 

0.364 

0.358 

1.338 

1.348 

1.334 


Clamped - Clamped 


0 

55.070 

55.189 

55.189 25.327 

25.336 

25.336 

90 

3.135 

3.145 

3.145 6.260 

6.264 

6.264 

[0/90], 

44.716 

44.805 

44.805 22.672 

22.679 

22.679 

[45/ - 45], 

5.478 

5.576 

5.495 8.275 

8.341 

8.280 


Simply- Supported 


0 

18.304 

18.307 

18.307 

13.430 

13.431 

13.431 

90 

0.812 

0.813 

0.813 

2.829 

2.829 

2.829 

[0/90], 

15.689 

15.692 

15.692 

12.434 

12.435 

12.435 

[45/ - 45], 

1.419 

1.496 

1.492 

3.739 

3.840 

3.834 


Exi/Eyy = 25, Gxy/Eyy = 0.5, Gy Z fEyy = 0.2, G XZ /Eyy - 0.5, V X y - 0.25, h 0 /b„ 0.5, 

Eyy = 1.9584 x 10 8 psf, bo = 0.1 ft, p = 0.250387 slugs/ft 3 , 

(E yy = 9.37687 x 10 9 Pa, b Q = 0.03048 m ,p = 129.0290 kg/m 3 ) 
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Table 4.4: Dimensionless fundamental frequencies (<2>) and dimensionless buckling loads 
(A) for symmetrically laminated beams with i/h a = 100 (a = 0, implies warping is not 


included; a = 1, imp 

.ies warping is included) 


A 

d) 

e 

Reddy 

Present 

Reddy 

Present 

(deg) 

(1997) 

(a = 0) 

(a = 1) 

(1997) 

(a = 0) 

(a = 1) 

Clamped-Free 






0 

5.134 

5.134 

5.134 

5.070 

5.070 

5.070 

90 

0.205 

0.206 

0.206 

1.015 

1.015 

1.015 

[0/90] s 

4.525 

4.525 

4.525 

4.758 

4.758 

4.758 

[45/ - 45] s 

0.358 

0.364 

0.359 

1.341 

1.352 

1.341 

Clamped - Clamped 






0 

80.665 

80.908 

80.908 

31.899 

31.916 

31.916 

90 

3.283 

3.294 

3.294 

6.450 

6.453 

6.453 

[0/90] s 

70.748 

70.969 

70.969 

29.857 

29.873 

29.873 

[45/ - 45] s 

5.737 

5.847 

5.755 

8.526 

8.599 

8.531 

Simply-Supported 






0 

20.461 

20.465 

20.465 

14.210 

14.211 

14.211 

90 

0.822 

0.822 

0.822 

2.848 

2.849 

2.849 

[0/90], 

18.015 

18.019 

18.019 

13.334 

13.335 

13.335 

[45/ - 45], 

1.436 

1.516 

1.510 

3.765 

3.867 

3.861 


E xx /E yy = 25, G xy /E yy = 0.5, G yz /E yy = 0.2, G xz /E yy = 0.5, u xy = 0.25, h o /b 0 = 0.5, 
E yy = 1.9584 x 10 8 psf, b 0 = 0.1 ft, p = 0.250387 slugs/ft 3 , 

(E yy = 9.37687 x 10 9 Pa, b Q = 0.03048 m, p = 129.0290 kg/m 3 ) 
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Table 4.5: 
with tjh 0 
included) 


Dimensionless fundamental frequencies (u>) for symmetrically laminated beams 
= 60 (a = 0, implies warping is not included; a = 1, implies warping is 


6 

Maiti k Present 

(deg) 

Sinha (1994) a = 0 a = 1 


Clamped-free 



[0/30/0] 

12.862 

12.869 

12.858 

[0/45/0] 

12.794 

12.796 

12.792 

[0/60/0] 

12.770 

12.769 

12.769 

[0/90/0] 

12.778 

12.778 

12.778 

Clamped - Clamped 



[0/30/0] 

80.753 

80.418 

80.354 

[0/45/0] 

80.282 

79.923 

79.902 

[0/60/0] 

80.065 

79.703 

79.703 

[0/90/0] 

80.047 

79.692 

79.692 

Simply- Supported 



[0/30/0] 

36.075 

36.051 

36.045 

[0/45/0] 

35.869 

35.840 

35.839 

[0/60/0] 

35.792 

35.760 

35.760 

[0/90/0] 

35.813 

35.780 

35.780 


E xx /E yy = 13.7088, G xy /E yy = 0.5471, G yz /E yy = 0.269641, G xz /E yy = 0.45679, 
v xy = 0.3, h 0 jb 0 = 0.3175, E yy = 9.42512 x 10 9 Pa, b D = 0.01 m, p = 1550.0666 kg/m 
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Table 4.6: Dimensionless fundamental frequencies (w) for unsymmetrically laminated 
beams with i/h 0 — 60 (a = 0, implies warping is not included; a = 1, implies warping is 


included) 


e 

(deg) 

Maiti & 
Sinha (1994) 

Present 
a = 0 a = 1 

Clamped- free 




[0/90/0/90] 

8.843 

8.853 

8.853 

[0/30/ - 30/0] 

12.397 

12.403 

12.403 

[0/45/ - 45/0] 

12.268 

12.270 

12.270 

Clamped- Clamped 



[0/90/0/90] 

56.000 

55.753 

55.753 

[0/30/ - 30/0] 

77.966 

77.580 

77.580 

[0/45/ - 45/0] 

77.086 

76.703 

76.703 

Simply- Supported 



[0/90/0/90] 

26.378 

26.352 

26.352 

[0/30/ - 30/0] 

34.786 

34.745 

34.745 

[0/45/ - 45/0] 

34.402 

34.370 

34.373 


E xx /E yy = 13.7088, Gxy/Eyy = 0.5471, G yz /E yy = 0.269641, G X z/E yy = 0.45679, 
i/ = 0.3, ho/b 0 = 0.3175, E yy = 9.42512 x 10 9 Pa, b a = 0.01 m, p = 1550.0666 kg/m 


4.5. SUMMARY 
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Table 4.7: Dimensionless buckling loads (A) for symmetrically laminated beams with 
(,/h 0 = 60 (a = 0, implies warping is not included; a — 1, implies warping is inclu ded) 

6 Clamped-free Clamped-Clamped Simply-Supported 

(deg) a = 0 a = 1 ct = 0 ot = 1 c* = 0 a = 1 

Unidirectional 

[0] 33.756 33.756 525.58 525.58 134.22 134.22 

[30] 10.026 5.633 158.96 90.343 47.158 39.989 

[45] 4.339 3.454 69.254 55.222 18.573 17.648 

[60] 2.850 2.736 45.566 43.750 11.514 11.424 

[90] 2.467 2.467 39.450 39.540 9.8610 9.8610 

Symmetric 

[0/0/0] 33.756 33.756 525.58 525.88 134.22 134.22 

[0/30/0] 33.120 33.063 515.33 515.33 131.71 131.67 

[0/45/0] 32.744 32.726 509.05 508.79 130.16 130.15 

[0/60/0] 32.607 32.607 506.44 506.44 129.59 129.59 

[0/90/0] I 32.652 32.652 | 506.60 506.60 [ 129.74 129.74 

E xx /E yy = 13.7088, G X y/E vy = 0.5471, G yz /E yy = 0.269641, G xz /E vy = 0.45679, 
v = 0.3, h 0 jb 0 = 0.3175, E yy = 9.42512 x 10 9 Pa, b a - 0.01 m, p = 1550.0666 kg/m 3 


Chapter 5 
A Probabilistic Approach 


In the past, the effects of uncertainties were recognized using traditional ap- 
proaches. These approaches simplify the problem by considering the uncertain 
parameters as deterministic and account for uncertainties using empirical safety 
factors. However, the conventional methods of design and analysis are not appro- 
priate for problems involving innovative design because the factors of safety are 
based on experience and there is no experience available for these problems. 

In fact, laminated composite structures have inherent uncertainties involved in 
the manufacturing process, and the end product may have significant variations 
in properties around the mean values. Thus the uncertainties in material and 
geometric properties should be taken as random in the analysis. In this context, 
the safety factors cannot be properly considered. Moreover, these safety factors do 
not provide any information on how the different parameters influence the overall 
behavior of the structure. 

In the previous chapters we have presented a finite element formulation to study 
the deterministic stability and vibrational response of laminated beams. However, 
this formulation does not include uncertainties. Here we intend to expand the deter- 
ministic formulation to account for uncertainties. Because probabilistic models can 
capture the influence of these uncertainties, the chapter begins with a description 
of the probability approach. Then we proceed to describe the random variables and 
their characteristics in this work. The rest of the chapter is dedicated to explain 
the three probabilistic theories developed here: probabilistic finite element method, 
sensitivity-based Monte Carlo simulation, and Monte Carlo simulation. 
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5.1 Uncertain Models 

Uncertainties in laminated composites exist because of material defects such as 
interlaminar voids, delamination, incorrect orientation, damaged fibers, and vari- 
ation in thickness. According to the nature and extent of uncertainty existing in 
laminated composites, different approaches can be used. 

If the uncertainty is due to imprecise information and/ or statistical data cannot 
be obtained, then the non- probabilistic approaches such as fuzzy sets can be used. 
These approaches have been studied by Elishakoff et al. (2001). On the other 
hand, if the uncertain parameters are treated as random variables with known 
(or assumed) probability distributions, then the theory of probability or random 
processes can be used. 

Probabilistic models can capture the influence of noncognitive sources of uncer- 
tainty because they are based on probability principles rather than on experience. 
These principles are mainly based on the following three axioms (Papoulis, 1991): 
(i) the probability of any single event occurring is greater or equal to zero; (ii) the 
probability of the universal set is one; (iii) the probability of the union of mutually 
exclusive events is equal to the sum of the probabilities. 

Several probabilistic methods have been used to analyze an uncertain unsym- 
metrically laminated beam by integrating uncertain aspects into the finite element 
modelling such as the perturbation technique using Taylor Series expansion and 
simulation methods (e.g., the Monte Carlo Simulation). 

Vinckenroy et al. (1995) presented a new technique to analyze these structures 
by combining the stochastic analysis and the finite element method in structural 
design. However they did not extend their work to dynamic problems. Stochastic 
methods were also studied by Haidar and Mahadevan (2000a). They applied the 
concepts to reliability analysis using the finite element method. 
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The probabilistic methods, as used here, are an extension of stochastic methods. 
The main difference is that stochastic finite element methods consider spatially 
correlated random variables and in probabilistic methods these random variables 
are not necessarily spatially correlated. Here we use the probabilistic finite element 
method with uncorrelated random variables. 


5.2 Random Variables 

5.2.1 Definition 

A random variable is defined as an uncertain parameter, for example, ply-angles, 
ply-thickness, length of the beam, width, etc. In the present work, the random 
variables are considered as independent and are denoted as 

r = {ri,r 2 ,...,r n } (5.1) 

where r^s are the different random variables. In the present work, the random 
variables for the laminated structures considered here are 

1. each ply’s orientation, 8i 

2. each ply’s axial Young’s modulus, E XXi 

Since these are independent random variables, the probability density function 
can be expressed as follows: 

n 

/(n,r 2 ,...,r n ) = JJ/i(r) 

1=1 


(5.2) 
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The probability density function (PDF) does not provide information on the prob- 
ability but only indicates the nature of the randomness. Among the used density 
functions in the analysis of structures are the Beta distribution, Normal or Gaussian 
distribution, Lognormal distribution, and Weibull distribution (Vinckenroy, 1994). 
From these, the most commonly used distribution is Gaussian. Thus the present 
analysis will assume that all random variables obey a normal distribution: 



where of and /Xj are the variance and the mean value of the i th random variable, 
respectively. 


5.2.2 Function of Multiple Random Variables 


In problems where uncertainties are considered, there exists no density function 
describing the random nature of the system. The information is limited to only the 
mean values of the random variables. In such cases, perturbation techniques are 
suggested, among other existing techniques (Ang and Tang, 1975; Schueller, 1997). 

In general the relationship between the random variables r, s and the matrix Y 
can be represented as a function of random variables, 

Y = y(ri,ra,...,r n ) (5-4) 


In most cases the sensitivity derivatives of matrix Y can be obtained. Therefore, 
the matrix Y can be expanded using Taylor series expansion about the mean values 
(Ang and Tang, 1984): 


Y(r u r 2 , . • ■ , r n ) = Y° + ^ Yjei + ~ + ’ ' 


j=i 


i=\ j = l 


(5.5) 
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&Y | 

dTidTj !x=x° 


where x° = (fii, /x 2 , . . . , /i„) is a set of mean random variables and ^ is a 

set of zero-mean uncorrelated random variables. 


5.2.3 Random Number Generator 


Most computers have the capability to generate uniformly distributed random num- 
bers Ui between 0 and 1. Afterwards random variables r, are, in general, obtained 

as 

s, = Si( ri ) = 4T 1 ( Ui ) (5-6) 

where s* is a function of the random variable r t , and < f >_1 is the inverse of the 
cumulative distribution function. However, when working with normal distributed 
random variables, one faces the challenge of not having the inverse of the normal 
distribution function in a simple closed-form expression (Law and Kelton, 2000). 
Various random generators exist, among these those developed by Odeh and Evans 
(1974) and Atkinson and Pearce (1976). In this dissertation, we use a slight mod- 
ification of the Box Muller transformation (Press et al., 1986). Basically, we start 
with two independent random numbers, u\ and u?, which come from a uniform dis- 
tribution (in the range from 0 to 1). Then we apply the Box Muller transformation 
to get two independent random variables which have a Gaussian distribution with 
zero mean and a standard deviation of one. Here we have used the modified version 
of the transformation which can generate random variables for any given mean and 
standard deviation. We have simplified the subroutine written in FORTRAN and 
is given as follows: 

SUBROUTINE randomvar (Nsamp .FirstSTD , Niuean , Randvar ) 

IMPLICIT none 
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! Parameters 

DOUBLEPRECISION, PARAMETER :: zero=0.0D+0, one=1.0D+0 

I 

! Used variables in this program 
INTEGER, INTENT (IN) :: Nsamp 

! 

! Random variables 

DOUBLEPRECISION, INTENT (IN) :: FirstSTD, Nmean 
DOUBLEPRECISION, INTENT (OUT) :: Randvar (Nsamp) 
DOUBLEPRECISION :: rand, xl ,x2,ul ,u2,Randvar2 (Nsamp) ,pi 

! 

! Counts random numbers 

INTEGER :: ii, jj, kk 

! 

! Holds the hour, minute, and second 
INTEGERS timeArray (3) 

! 

pi = dacos(-one) 

i j t !! I ! I j j 1 ! M m i ! j j !!!!!! m !!!!!!!!!!!!!!!!!!!!!!!!!!!! ! 
! ! 

! Initializing the seed of the random number function ! 
! with the sum of the current hour, minute, and \ 

! second to get a different sequence each time \ 

I ■ 

M I I j M m m M M i M j | I m jj!!!!!!!!!!!!!!!!!!!!!!!!!! M !!! ! 

call it ime (timeArray) ! Get the current time 

ii = rand( timeArray(l)+timeArray(2)+timeArray(3) ) 

j j j j j j i j j m j j mi j j j j j ; j j !!!!!!!!!!!!!!! ! ! ! ! ! ! ! !!!!!! M ! ! 

j 

j Calling rand(ii) , the numbers are generated 
! between 0 and 1 . 

! 

! log(u) -> returns the natural logarithm of u 

j 

j j j j j j ; j j j j j m j j j j j j i ; ; ; j j (jj!!! m !!!!!!!!!!!!!!!!!!!! ! ! 
Do kk=l, (Nsamp/2) 

ul = rand(ii) 
u2 = rand(ii) 
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xl = FirstSTD*dsqrt( -2.0D0*dlog(ul) )*dcos(2.0D0*pi*u2) 
& + Nmean 

x2 = FirstSTD*dsqrt( -2.0D0*dlog(ul) )*dsin(2.0D0*pi*u2) 
& + Nmean 

Randvar ( kk ) = xl 

Randvar ( (Nsamp/2)+ kk ) = x2 

EndDo 

return 

END SUBROUTINE randomvar 


5.3 Monte Carlo Simulation 

Monte Carlo Simulation, although computationally expensive, is a quite versatile 
technique that is capable of handling situations where other methods fail. The 
MCS is also often used to verify the results obtained from other methods. 

Monte Carlo Simulation methods are based on the use of random variables and 
probability statistics to investigate problems. Vinckenroy and de Wilde (1995) 
developed an approach to account for uncertanties by combining the finite element 
method and Monte Carlo techniques. 

A large sample is generated and then using probability density functions (PDF) 
one evaluates the probability of having such values. The larger the number of 
simulations, the higher the confidence in the probability distribution of the obtained 
results. Therefore, for the present analysis, at least ten thousand realizations of the 
uncertain beam are performed, increasing the accuracy of the ply-angle and ply 
axial modulus of elasticity distribution fit to the sample data. 
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5.4 Probabilistic Finite Element Analysis 


Let K be the probabilistic dimensionless linear stiffness matrix, M the probabilis- 
tic dimensionless mass matrix, L the probabilistic dimensionless loading matrix, A*, 
the probabilistic dimensionless eigenfrequency of the k th mode, Pk the probabilistic 
dimensionless buckling load of the fc th mode, rj) k the probabilistic dimensionless left 
eigenvector of the k th mode, and 4> k the probabilistic dimensionless right eigenvec- 
tor of the k th mode. Thus the probabilistic dimensionless eigenvalue problem is 
expressed as 

{r!> k } T [K-P k L-\ k M} = 0 (5J) 

[K-P k L-X k M]{ ( j> k } = 0 
Note that the above quantities are all dimensionless. 

In the stability and free vibration analysis of the present problem, the random 
nature of the stiffness matrix, mass matrix, loading matrix, buckling loads, eigenfre- 
quencies, and eigenvectors are studied using a Taylor series expansion up to second 
order about the mean of each random variable. In contrast to the work done by Oh 
and Librescu (1997), the present formulation can use results provided by commer- 
cial finite element codes, greatly reducing the number of calculations. The present 
formulation also calculates the eigenfrequency sensitivities up to second-order ap- 
proximation. Moreover, the method is extended for nonconservative systems and 
for the analysis of stability of laminated composites using the dynamic criterion. 

The presence of structural uncertainties results in a certain randomness in the 
extensional matrix A, extensional-bending coupling matrix B, and bending matrix 
D. The coefficients of these matrices are present in the equivalent bending stiffness 
matrix D c , Eq. (2.37). Thus the matrix D c will have certain randomness associated 
as well and, by virtue of Eq. (3.38), these uncertainties affect the stiffness matrix. 
Thus, the random nature of the stiffness matrix K is expanded in terms of the 
mean-centered zeroth-, first-, and second-order rates of change with respect to the 
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random variables, as presented in Eq. (5.5), 


K(r u r 2 , . . . ,r n ) = K° + K 1 ^ + ^ E K ” e * e i 


(5.8) 


i=l 


i = 1 j=l 


Similarly, the random natures of the mass matrix M and the loading matrix L are 
expanded in terms of their mean-centered zeroth-, first-, and second-order rates of 
change with respect to the random variables as 


n 1 n n 

Af(ri,r 2 ,.. . ,r n ) = M° + ^ M'e, + 5 EE M « 

1=1 i=l i=l 

l[t 1, 7*2, , r" n ) = i' 0 + y: + g y: 


i=l 


1 = 1 j= 1 


(5.9) 

(5.10) 


The buckling loads, eigenfrequencies, and left and right eigenvectors are also af- 
fected by uncertainties. The eigenfrequencies and eigenvectors are expressed in 
terms of their mean-centered zeroth-, first-, and second-order rates of change with 
respect to the random variables. Further, let the first-order rate of change of the 
k t ' 1 mode right eigenvector with respect to the i 1 * 1 random variable evaluated about 


the mean be 


dx i 



(5.11) 


and the second order rate of change of the k th mode right eigenvector with respect 
to the z th and j th random variables evaluated about the mean be 


dxi dxj 



(5.12) 


Using the same notation for the eigenfrequencies and buckling loads, we can express 
the following: 
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Eigen frequencies 


n ^ n n 

A*(ri, r 2 , • • • , r n ) = \° k + ^ Ajy e { + - ^ ^ 


i=l 


t=l J — l 


(5.13) 


Right Eigenvectors 


{<t>k( r ur 2 , . ■ 


r n)} = {<Pk} + + 

1=1 




(5.14) 


Buckling Loads 


Pk(r u r2,. 


■ .r„)=P° t + 'EPl,^ + l'E't, P ^ 


II ff 

kij e » e J 


(5.15) 


t= 1 


i= 1 j=l 


The substitution of Eqs. (5.8), (5.9), (5.10), (5.13), (5.15), and (5.14) into 
Eq. (5.7), results in a probabilistic eigenvalue problem. Since the uncertainties 
in the random variables are assumed small, in the applied perturbation technique 
it is sufficient to only consider up to second-order terms. Thus the expansion of the 
probabilistic eigenvalue problem leads to three equations which are solved succes- 
sively: 


/’ : [K° - P° t L° - \tM«} {<*.“} = 0 

(5.16) 

e‘ : [Kl - P' ki L° - (<*£} = 0 

(5.17) 

£ 2 : [Kij - P" £° - A&M 0 ] M = 

(5.18) 

- [JT' - P' U L° - A^M 0 ] «} 



Details are included in Appendix E. 
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5.4.1 Eigenfrequency Derivatives 


Equating the zeroth order terms of e* in the eigenvalue expansion, an eigenvalue 
problem for the mean- valued system is obtained. Therefore, the mean-centered ze- 
roth derivative eigenfrequencies and associated eigenvectors are obtained as follows: 


[K° - P° k L° ~ A° fc M°] = 0 (519) 

[K° - P° k L° - A°M°] {$} = 0 

Recall that the loading and mass matrices do not depend on the lamina s mechan- 
ical characteristics. Since the lamina ply angles and axial Young s modulus are the 
only random variables considered, the sensitivity derivatives of the mass and load- 
ing matrix vanish. Thus expressions for the mean-centered first and second-order 
derivatives of the eigenfrequencies are found as 


[k[ - pj£] {#} 

* " {<} T M 

„ [Kj[ - P&£°] {&} 

*« {*J} T [M 1 ] {<( 

(V’£) T {K[ - PIT 1 - A^M"] «} 

+ {WWTM 

+ M T [ m °] M 


For conservative systems, 


L = L t and ° k } = {< Pi } 

Thus for conservative systems, and by virtue of Eq. (5.17), the expression for the 
second-order derivative simplifies to the expression derived by Kapania and Goyal 
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(2002) for conservative systems: 

xI1 _ [Kl[ - P"L°] { 0 ° fc } 2 

{<t>lY [M°] M 

The advantage of this method is that the eigenvalue problem needs to be solved 
only once. The sensitivity analysis is done by using results from the mean- valued 
eigenvalue problem. This results in great computational saving. 


5.4.2 Eigenvector Derivatives 

Although in this dissertation we have not studied the derivatives of eigenvectors, 
the method on how to determine these derivatives is described in Appendix E, and 
in this chapter we only highlight the final derivation. Recall that to normalize 
the left and right eigenvectors we used two independent criteria. Since the right 
eigenvectors form a complete set of vectors, an eigenvector can be represented by 
a linear combination of all other eigenvectors: 

{*9 = IMS {*4 < 5 ' 23 > 

3 = 1 


where 


a {i) - 
a kj ~ 


A*, — A j 


j = k 


j 7 ^ 


(5.24) 
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5.4.3 Stiffness Matrix Derivatives 

Hasselman and Hart (1972) proposed a method to calculate derivatives of eigen- 
frequencies for reduced systems. However, they assumed that the contribution of 
the partial derivatives of the transformation matrix to the partial derivatives of the 
eigenfrequencies and eigenvectors is small. 

In calculating derivatives of eigenfrequencies, only derivatives of the stiffness and 
mass matrices are needed. The derivatives of the stiffness matrix are more involved 
because they require taking the derivatives of the equivalent bending-stiffness ma- 
trix, Eq. (2.37). Various numerical schemes, such as the finite difference method, 
exist to evaluate these derivatives. When using some of these numerical schemes, 
ill-conditioning could be a problem. This problem can be avoided by the follow- 
ing formulation which allows the derivatives to be obtained exactly by numerical 
multiplication. The technique consists of taking the derivatives of Eq. (2.36) with 
respect to each of the variables xj 

[D R \ Ti = [£>/,/] >r . - [Diiji ]' 1 [DjiA (5-25) 

— [Diji] [D 11 ,11] A. [Dii,i] ~ [Diji] [D n,n] [Diij] r . 


[D R ], ri , 


[DjA, rirj ~ [A;,//], rir . [Diiji] 1 [Du A (5-26) 

- [Dull, [Du, nl) ~ l D ^u\~ l [■ Diul , 


- \Duil, [Diiji] r ) [D n j] ~ [Di,ii] [D u,ii\, T ) rj [Dii,i] 

- [Diji] [Diiji ]' 1 [Dnj] r . - [Diji] r . [Diiji] 1 [Diij] r . 

- [Diji] [Diiji] r ) [Dnj] r . — [ Diji ] [Diiji] [Dnj] r . r . 


The derivatives of [ Djjjj ] 1 are calculated using the following matrix definition. 

[Diiji ]' 1 [Diiji] = [/] (5-27) 
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and these derivatives are 


{DllJiti = 

- [Dii,ii] 1 [Dn,u) r . [Dnji] 1 

(5.28) 

[Diiji]']., = 

- [Pnji]-) [ £>//,//], r 4 [^//, z/] -1 

(5.29) 


- [£>//,//] _1 [-0//,//], r , [■ D //.//], r , 1 



- \Dn, //] _1 [- D //,//],r i r j 



5.5 Statistical Quantities 


The mean value of the eigenfrequency is obtained by taking the expected value of 
Eq. (5.13): . „ 

» = EM = A2 + i £ E E ^ (5 ' 30) 

i=l j = 1 

The variance and the standard deviation are defined as 

Var[A fc ] = E[X 2 k ] - (5-31) 

<7a = \ZV^ar[A*;] (5.32) 

The uncertainty associated with the inherent randomness is given by the coefficient 
of variation (Ang and Tang, 1975) 


6\ = a\/v\ 


(5.33) 


For symmetrically distributed independent random variables, 


Var[\ k ] = £>LaL £[<?] + \±± A2.A" {E[e?€?] - «[«} (5.34) 


i= 1 


i=l jr = 1 


5.6. SUMMARY 


147 


where 

m 

Em 

E\A\E[<*\ 

and Ngam is the number of samples. 

5.6 Summary 

The sources of uncertainties used throughout this work are those for which un- 
certain parameters can be treated as random variables with known (or assumed) 
probability distributions. Here we assume all random variables to be independent 
and spatially uncorrelated. In this chapter, we develop the probabilistic finite ele- 
ment model, which is based on the deterministic response. The eigenvalue problem 
has to be solved only once, this being a great computational advantage. Sensitivity 
derivatives for conservative and nonconservative systems are presented as follows: 

At = {<) T [*' - -Pit 0 ] «} ( 5 - 35 > 

A& = {<} T m - P&L°] {<*>2} (5.36) 

+ {^2 f [Ki - P' U L° - a2,M°] {<} 

+ W2r [*J - p ki L ° - A l,M°] {<*>L> 


Aaam / \2 

E V g ~ l*i) 

N sam 

9=1 

( r 9 ~ ) 2 ( r g ~ Mj) 2 

N 

. x y sam 

9=1 

Asam f* / , \ / „ , . \ 

V"^ V""' \ r q ~~ f^i ) v r ~~ MjJ 

, . ^sam 

o=l r=l L 


In the next chapter we perform the analysis of uncertain laminated beams. 


Chapter 6 

Reliability Analysis of Laminated 

Beams 


The analysis in Chapter 4 is only valid for perfect structures, those structures 
without imperfections. However, uncertainties in ply angles and the axial modulus 
of elasticity lead to imperfections in the structure and the deterministic analysis 
may be no longer valid. Thus here we use the methods described in the Chapter 5 
to perform a more accurate analysis. The analysis of uncertain laminated beams is 
performed to two types of problems: free vibration analysis and stability analysis 
using the dynamic criterion. 


6.1 Probabilistic Analysis 

In a problem involving uncertainty, one first conducts statistical analysis on the 
random variables. This can be obtained experimentally or using sampling tech- 
niques. Then using this information one calculates the influence of the randomness 
of the random variables on the wanted response. 

In the present analysis, we used ten thousand data points and the results are 
presented in frequency density diagrams or histograms, which show the distribu- 
tion of the eigenfrequencies or buckling loads. The number of cells needed in the 
frequency density diagram is given by Sturges’ formula as (Law and Kelton, 2000): 
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Jfcfd « 1 + 3.3 log 10 (iV 8am p) (61) 

where iV samp is the number of data points. For N aamp = 10, 000, the approximate 
number of cells is 14.2. However, we have chosen to use fcfd = 16 in this work. Once 
we have the histograms, we pick a density function that best fits the response. 
This probabilistic density function can be used to perform the reliability analysis 
of uncertain structures. 


6.1.1 Random Variables 


In the present study, it is assumed that the beam is composed of identical plies that 
possess the same geometric and mechanical properties. It is further assumed that 
the randomness of each ply angle and modulus of elasticity in the ^-direction, E XXi , 
is the same for every ply and are spatially uncorrelated. Let 6 1 and E XXt be the 
deterministic quantities of the i th lamina. If the probability distribution function 
for the ply angles and the axial modulus of elasticity are f e (r) and f E {r), then the 
randomness can be expressed as 


@i = + 0r 

E XXi = &xxi H" Er 

where i is represents the i th ply, and 0 T and E r represent the random variables 
expressing the uncertainty in §i and E XXi , respectively. 

In general, the ply-angle uncertainties are between ±2.5°. Thus for the present 
study we have assumed for 6 r a Gaussian distribution with a standard deviation 
of 2.5°. Thus there is a 95.3% probability that the ply orientation will have an 
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0.18 



Ply Angle, 6 (deg] 


Figure 6.1: Probability density function for the full generation of possible variation of the 
ply angle, 6 , with a standard deviation of oq — 2.5° and zero mean. 


uncertainty between —5° and 5°: 

P( -5° < r < 5°) = / 1 7 = r * ) dr = 0.953 (6.2) 

J-v *0 \FT^ 

The probability density function is shown in Fig. 6.1. 

The randomness in the material properties, with a 95% confidence interval, have 
an experimental coefficient of variation of 3% (Lin and Kam, 2000). However, in 
the present work we have assumed a coefficient of variation of 5%. Thus there is 
a 95% probability that the ply orientation will have an uncertainty between -0.1 
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. 0,19 -0.16 -0.13 -0.1 / -0.08 -0.05 - 0.03 0.00 0.03 0.05 0.08 Q.ll 0.13 0.16 0.18 0.2 i 


Dimensionless Axial Modulus of Elasticity, n xv 


Figure 6.2: Probability density function for the full generation of possible variation of the 
dimensionless axial modulus of elasticity for isotropic: beams, n ss — 1, with a dimension- 
less standard deviation of a f? = 0.05 and zero mean. 
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0.60 



Dimensionless Axial Modulus of Elasticity, n ^ 


Figure 6.3: Probability density function for the full generation of possible variation of 
the dimensionless axial modulus of elasticity for laminated beams, n xs — 13.7088, with a 
dimensionless standard deviation of rr# = 0.69 and zero mean. 
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and 0.1: 

P(-0.1 < r < 0.1) = [ 1— ; e ~^ dr = 0.953 (6.3) 

v J- o.i a E s/Tir 

The probability density function for the various cases studied in this dissertation 
are shown by Figs. 6.2 and 6.3. Note that because we have nondimensionalized all 
quantities in the finite element formulation, we take the variation of the nondimen- 
sional axial modulus of elasticity: 


6.1.2 Probabilistic Models 

For the uncertain analysis of laminated beams three models are developed. Ex- 
act Monte Carlo Simulation (EMCS), Sensitivity- Based Monte Carlo Simulation 
(SBMCS), and Probabilistic Finite Element Analysis (PFEA). The analysis can be 
described as follows: 

1. In the EMCS we solve the eigenvalue problem for each realization to determine 
the natural frequencies and buckling loads. The mean value and the standard 
deviation are obtained using statistical methods. 

2. In the SBMCS, we only solve the eigenvalue problem for the mean values of 
the random variables. Then using this information, we calculate the sensi- 
tivity derivatives, leading to an approximate function describing the random 
nature of the eigenfrequencies and buckling loads. Now we calculate the 
eigenfrequency and buckling load for each realization. The mean value and 
the standard deviation are obtained using statistical methods. 

3. The PFEA uses the sensitivity derivatives and gives the mean and standard 
deviation directly. 
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6.1.3 Finite Element Analysis 

For all three methods above described we used the finite element method with 
five finite elements. The three sets of boundary conditions used are the same as 
those mentioned in Chapter 4: hinged-hinged, clamped-free, and clamped-clamped. 
Moreover, recall that all the analysis was performed using the dimensionless finite 
element method, and the response is in nondimensional form. When studying the 
effect of uncertainties of random variables on the fundamental natural frequencies, 
it is convenient to study their squared value, i.e., eigenfrequencies, which are given 
in their nondimensional form as 

(6.4) 


(6.5) 

6.2 Free Vibrations 


An — A n 


Jo l 4 

Eyy b Q h Q 


A ji A n 


12 ip l 4 

Eyy bo h Q 


The critical loads are given as 


P = P 

1 n — 1 n 


Eyy b 0 h Q 


Pn = P n 


12 r 

Eyy bo h a 


We first study the influence of having an uncertain dimensionless Young’s modulus, 
n xx , on the free vibrations of isotropic beams. Next, we study how the free vibration 
response of various laminated beams are affected by uncertainties in n xx . Also, we 
studied the cases for ply-angle variations. 
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For the case of free vibrations the sensitivity derivatives of the natural eigenfre- 
quencies simplify to: 


At = m T [*,'] 

(6.6) 

ai' = m T [k!‘] M) 

(6.7) 


6.2.1 Isotropic Beams: Uncertain Young’s modulus 


The probability distribution functions of the fundamental dimensionless eigenfre- 
quencies, A, are shown in Figs. 6.4, 6.5, and 6.6. Results show that by randomly 
generating possible values for the Young’s modulus in the ^-direction with a coeffi- 
cient of variation of 5%, the fundamental dimensionless eigenfrequencies also have a 
coefficient of variation of 5%. The figures show that a symmetric randomness in the 
random variables produces symmetric variation in the dimensionless fundamental 
natural frequency. Also, the Sensitivity-Based Monte Carlo Simulation (SBMCS) 
when using only one thousand samples are in perfect agreement to those by the 
Exact Monte Carlo Simulation method (EMCS) using ten thousand samples. 

It seems that the variation in E xx is most critical for a fixed-fixed isotropic 
beam. This could be because the boundary conditions make the beam stiffer and 
thus increases the fundamental frequency. Although the coefficient of variation is 
only of 5%, it significantly affects the fundamental frequency because of the high 
frequencies. 

For the case of Sensitivity-based Monte Carlo Simulation, we also studied the 
influence of the zeroth, first, and second order variation on the dimensionless natural 
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eigenfrequencies, where these orders are understood as follows: 


A* 



zeroth order 


+ Afc, e.e, 


i=l j — 1 


second order 


( 6 . 8 ) 


Figures 6.7, 6.8, and 6.9 show that the second order terms have no significant 
influence on the overall dimensionless fundamental eigenfrequency. 


6.2.2 Laminated Beams: Uncertain Young’s modulus 


Results for a unidirectional laminated beam with a ply angle of 0° are plotted in 
Figs. 6.10 and 6.11. Similar trends to those of the isotropic case can be observed. 
However, for a 90° unidirectional laminated beam, Figure 6.12 shows that the 
Young’s modulus in the x-direction has very little effect on the free vibrational 
response. The reason for this could be that the laminate in stiffer in the x-direction, 
thus small variations in E xx do not influence the beam s fundamental frequency. 
Results for other boundary conditions were consistent with those of the cantilevered 
case. 

We also studied several laminated composites such as sandwiched laminas with 
a layout of [O°/0/0/O°], for a11 6 = 30°, 90°. Figures 6.13-6.18 show these results. 
For all three boundary conditions it can be seen that the variation in the dimen- 
sionless fundamental frequency is the smallest for 0 = 90°. Once again a very good 
agreement holds for all three methods employed in this study. 



Probability density function 
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Dimensionless Eigenfrequency 


(a) Monte Carlo Simulation 



Dimensionless Eigenfrequency 


(b) All three methods used here 

Figure 6.4: Probability density function of the dimensionless eigenfrequency, A, for a 
simply-supported isotropic beam with uncertain Young’s modulus in the x-direction. 
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Dimensionless Eigenfrequcncy 


(a) Monte Carlo Simulat ion 



9.50 10.50 tl.50 12.50 1150 U50 

Dimensionless Eigenfrequency 

(b) All three methods used here 


Figure 6.5: Probability density function of the dimensionless eigenfrequency, A, for a 
cantilevered isotropic beam with uncertain Young’s modulus in the x-dircctioii. 
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(a) Monte Carlo Simulation 



Dimensionless Eigenfrcquency 

(b) All three methods used here 

Figure 6.6: Probability density function of the dimensionless eigenfrcquency, A, for a 
fixed-fixed isotropic beam with uncertain Young’s modulus in the ./’-direction. 
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Figure 6.7: Effect of the order of the sensitivity derivatives on the dimensionless eigenfre- 
quency, A, for a simply-supported isotropic beam with unc ertain Young’s modulus in the 
x-direction. 
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Figure 6.8: Effect of the order of the sensitivity derivatives on the dimensionless eigen- 
frequency, A, for a cantilevered isotropic beam with uncertain Young’s modulus in the 
x-direction. 
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Figure 6.9: Effect of the order of the sensitivity derivatives on the dimensionless eigen- 
frequency, A, for a fixed- fixed isotropic beam with uncertain Young’s modulus in the 
^-direction. 
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Figure 6.10: Probability density function of the dimensionless eigenfrequency, A, for a 
cantilevered unidirectional laminated beam (0°) with uncertain Young’s modulus in the 
x-direction. 
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Figure 6.11: Effect of the order of the sensitivity derivat ives on the dimensionless eigenfre- 
quency, A, for a cantilevered unidirectional laminated beam (0°) with uncertain Young’s 
modulus in the x-direction. 
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Figure 6.12: Probability density function of the dimensionless eigenfrequency, A, for a 
cantilevered unidirectional laminated beam (90°) with uncertain Young's modulus in the 
rr-direction. 
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Figure 6.13: Probability density function of the dimensionless eigenfrequency, A, for a 
fixed-fixed laminated beam ([0°/30 o /30 o /0°]) with uncertain Young’s modulus in the x- 
direction. 
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Figure 6.14: Probability density function of the dimensionless eigenfrequency. A, for a 
fixed-fixed laminated beam ( [0 o /90 o /90° /0°] ) with uncertain Young’s modulus in the x- 
direction. 
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Figure 6.15: Probability density function of the dimensionless eigenfrequency, A, for a 
cantilevered laminated beam ( [0° /30°/30 o /0°] ) with uncertain Young’s modulus in the 
^-direction. 
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Figure 6.16: Probability density function of the dimensionless eigenfrequency, A, for a 
cantilevered laminated beam ([0°/90 o /90 o /0°]) with uncertain Youngs modulus in the 
^-direction. 
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Figure 6.17: Probability density function of the dimensionless eigenfrequency, A, for a 
simply-supported laminated beam ([0°/30 o /30 o /()°]) with uncertain Young’s modulus in 
the :r-direction. 
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Figure 6.18: Probability density function of the dimensionless eigenfrequenry. A, for a 
simply-supported laminated beam ([0°/90 o /9() o /()°|) with uncertain Young’s modulus in 
the ^-direction. 
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6.2.3 Laminated Beams: Uncertain Ply Angles 

We also considered the cases when the ply orientations may become uncertain. For 
this case we studied several laminated composites such as sandwiched laminas with 
a layout of [O°/0/0/O°], for all 6 = 30°, 90°. Figures 6.19-6.22 show these results. 

In all three models, the mean values and the coefficient of variations were close. 
However, the probabilistic finite element analysis and the Sensitivity Based-Monte 
Carlo Simulation are conservative in the sense that both of them overestimate the 
variation of the natural frequencies. Exact Monte Carlo Simulations would have 
been the most accurate approach but also a very expensive one. Therefore, the 
probabilistic finite element analysis can be safely used. The Sensitivity-Based MCS 
is an alternative approach that produces fairly good results and saves time. This 
approach produces very good results for only one thousand samples as opposed to 
ten thousand samples employed in the exact MCS. 

It should be noted that the results presented by Kapania and Goyal (2002) had 
a better correlation among the three models. The reason is that in their work they 
considered ply variations between -5° and 5° only. Here we have decided to use 
the full spectrum of the data. However, in both cases it can be shown that the 
ply angle uncertainties can play an important role in affecting free vibrations of 
symmetrically and unsymmetrically laminated beams. 


6.3 Reliability Analysis 


In any given problem the variations may be highly significant or they may be 
insignificant. It is important to have an idea of their magnitudes. Sometimes this 
information is unknown and we use safety factors to overcome our lack of knowledge. 
In the sense of structural stability, a structure is safe only if the actual load applied 
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Dimensionless Eigenfrequcncy 

Figure 6.20: Probability density function of the dimensionless eigenfrequeney, A, for a 
cantilevered laminated beam ([() o /9() o /!)() o /0°|) with uncertain ply orientations. 
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Figure 6.22: Probability density function of the dimensionless eigenfrequency, A, for a 
simply-supported laminated beam ([0 0 /9() u /9() 0 /() 0 ]) with uncertain ply orientations. 
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to the component does not exceed the critical load. In the traditional method, the 
degree of safety is usually expressed by the safety factor, S F . and is defined as 

SF=^f (6.9) 

The factor of safety is used with the following criterion in mind: if SF = 1 then 
the structure is on the point of failure; if SF < 1 then the structure is in a failed 
state; and if SF > 1 then the structure is considered safe. 


A higher value of the safety factor seems to indicate a safer component. However 
this is not necessarily the case as the inevitable variations must be kept in mind. 
Let us consider the case of a structure having a critical load of and the nominal 
values of P = Pcr/SF. Now because of the inherent imperfection in the structure, 
let the probability density functions of the dimensionless critical load be /(A), as 
shown in Figure 6.23. It is our goal to find the probability of failure for a structure 
designed for various safety factors, i.e. , Ai = A cr/SF. 


In the present work, we calculate the normal distribution functions by using 
either the Exact MCS, sensitivity-based MCS, or the probabilistic FEA. It should 
be noted that we have assumed that the buckling loads obey a normal distribution, 
although in some cases this might not be true. Now we study the probability of the 
structure to fail under a given factor of safety. Let the probability density function 
have a mean /x p and a standard deviation o p . The probability of failure is then 
evaluated as follows: 


Vi 


-r 


/(A) d\ 


(6.10) 


The reliability of the system is defined as 1 — p/. 
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6.3.1 Conservative Case 


We first study the influence of having an uncertain dimensionless Young’s modulus, 
n xx , on the buckling of isotropic beams. Next, we study how the stability of various 
laminated beams is affected by uncertainties in n xx and ply angles. For the case of 
conservative buckling analysis, the critical load occurs when the eigenfrequency is 
zero. Thus the sensitivity derivatives of the buckling load simplify to 


„ 1M1 

* [L°] {</>° fc } 

l/; (0° fc } r [*ff] M 

kiJ [L°] {$} 


(6.11) 


( 6 . 12 ) 


For the case of isotropic beams, it was found that the dimensionless Young’s 
modulus in the x-direction does affect the dimensionless buckling load, although 
the variation is small. The variation of the dimensionless buckling load is shown in 
Figs. 6.24, 6.26, and 6.27. For the case of a fixed-fixed isotropic beam, the variations 
in the buckling load are most critical. 

We also studied various cases of unidirectional cantilevered laminated beams 
with a ply angle of 0°. Figures 6.28 and 6.29 show that ply angle variations signif- 
icantly affect unidirectional laminated beams when compared to variations in the 
dimensionless Young’s modulus in the x-direction. Also, it was found that E xx had 
no influence whatsoever on the buckling load of unidirectional laminated beams 
with a ply angle of 90°. A similar trend was found for all other unidirectional 
laminated beams. 

The probability of failure for a cantilevered isotropic beam with uncertain 
Young’s modulus in the x-direction is shown in Fig. 6.25. The reliability of the 
structure is shown as a second figure to all the cases mentioned above. 
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When performing the reliability analysis for all the conservative cases studied 
here, results showed that, for the uncertainties considered here, the structure is 
reliable when it is designed for a safety factor of 1.5, a value traditionally used 
in aerospace design. Thus, when uncertainties in ply orientations and Young’s 
modulus in the :r-direction affect the laminated beams, the structure can be safely 
modeled using deterministic approaches. 


6.3.2 Nonconservative Case 


For the nonconservative problem, we find the critical load such that the first two 
eigenfrequencies coalesce, known as flutter point. At the onset of flutter we know 


that 


dP 

d\k 


= 0 


(6.13) 


Recall that the zeroth order eigenvalue problem, given by Eq. (5.16), is defined as 


[K° - P° k L° - A °M°] {$} = 0 (6.14) 


Now premultiplying the above equation by the left eigenvector and solving for the 
critical load, we get 


i0 [K”] {$} Ml 
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(6.15) 


By virtue of Eq. (6.13), at the onset of flutter 


MY [M°] {$} = 0 


(6.16) 
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Thus the equations for the sensitivity derivatives for the critical load at the onset 
of flutter become 


“ {<) T [i°] {<^2} 


(6.17) 


wiv [Kin m (6 18) 

ki ‘ {^} T [£°1 {$} 

where the left and right eigenvectors, and the stiffness derivatives, are evaluated at 
the onset of flutter using the mean values of the random variables. 

Figure 6.30 shows the variation in the critical load, which occurs at flutter, 
for a purely tangential follower load. Results show that isotropic beams under 
nonconservative loading also have a small probability of failure. For the laminated 
beams, shown in Figs. 6.31 and 6.32, the structure will be reliable when it is designed 
for a safety factor of 1.5, a value traditionally used in aerospace design. 

Thus, as was the case for conservative loading, the laminated beams can be 
modeled using deterministic approaches for the type of uncertainties considered 


here. 
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(a) Probability density funct ion 



(b) Structure’s reliability against various safety factors 


Figure 6 . 24 : Probability density function of the dimensionless buckling load, P, and the 
structure’s reliability for a cantilevered isotropic beam with uncertain Young’s modulus 
in the ^-direction under a conservative compressive loading. 


Safety Factor 


6. 3. RELIA DILIT Y ANAL YSIS 


133 



Figure 6.25: Probability of failure cantilevered isotropic beam with uncertain Young's 
modulus in the x-direction under a conservative compressive loading. 
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Dimensionless Buckling Load 


(a) Probability density funct ion 



(b) Structure’s reliability against various safety factors 


Figure 6.26: Probability density function of the dimensionless buckling load, P, and the 
structure’s reliability for a fixed-fixed isotropic beam with uncertain Young’s modulus in 
the ^'-direction under a conservative compressive loading. 
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(a) Probability density function 



(b) Structure’s reliability against various safety factors 

Figure 6.27: Probability density function of the dimensionless buckling load, P , and 
the structure’s reliability for a simply-supported isotropic beam with uncertain Young’s 
modulus in the x-direction under a conservative compressive loading. 
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(a) Probability density function 



Reliability (%) 

(h) Structure’s reliability against various safety factors 


Figure 6.28: Probability density function of the dimensionless buckling load, P, and the 
structure’s reliability for a cantilevered unidirectionally laminated beam of a ply of 0° and 
uncertain Young’s modulus in the .r-direction under a conservative compressive loading. 
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(b) Structure’s reliability against various safety factors 

Figure 6.29: Probability density function of the dimensionless buckling load, P, and the 
structure’s reliability for a cantilevered unidirectionally laminated beam of a ply of 0° 
and uncertain ply angle under a conservative compressive loading. 
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Dimensionless Buckling Load 


(a) Probability density funct ion 



(b) Structure's reliability against various safety factors 


Figure 6 . 30 : Probability density function of the dimensionless buckling load, P, and the 
structure’s reliability for an isotropic cantilevered beam with uncertain Young’s modulus 
under a purely tangential follower load. 
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Dimensionless Buckling Loud 

(a) Probability density function 



(b) Structure’s reliability against various safety factors 


Figure 6.31: Probability density function of the dimensionless buckling load, P, and 
the structure’s reliability for an undirectional laminated cantilevered beam (0°) with 
uncertain Young’s modulus under a purely tangential follower load. 
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(u) Probability density function 
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(b) Structure’s reliability against various safety factors 


Figure 6.32: Probability density function of the dimensionless buckling load, P, and 
the structure’s reliability for an undirectional laminated cantilevered beam (0°) with 
uncertain ply angle under a purely tangential follower load. 


Chapter 7 
Concluding Remarks 


Fiber-reinforced laminated composites have found widespread use in aerospace 
engineering. Moreover, these composites can exhibit a considerable variation in 
their material properties because their manufacturing process involves a number of 
parameters that may not be fully controllable. 

In the aerospace industry, the stability analysis of laminated composite struc- 
tures is extremely important. In this work, we have studied the stability of lami- 
nated structures subject to subtangential loading using the dynamic criterion. We 
also studied the influence of uncertain ply orientations and uncertain axial Young s 
modulus on the stability of these laminated structures. 

In this chapter, we intend to present our final remarks for the deterministic 
analysis as well as the probabilistic one. The last section is devoted to present 
various topics to which this work can be extended. 


7.1 Deterministic Response 

A twenty-one degree of freedom laminated beam element has been developed for the 
one-dimensional analysis of symmetrically and unsymmetrically laminated compos- 
ites. The element has the following features: use of symbolic integration, all three 
displacements (axial, transverse, and lateral), torsional and warping effects, inplane 
rotation, and shear rotation. Here we used the simplest possible warping function, 
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one that would satisfy the classical laminate theory. The entire derivation is pre- 
sented in its dimensionless form, which often results in a computationally robust 
approach. An advantage of this method is that the coefficients of the stiffness 
matrix are of the same order. The same is true for the mass matrix. 

The element takes into account the various couplings that affect laminated com- 
posites, such as inplane shear and twisting rotation. Although the inertia due to 
shear rotation is usually small, its relevant contribution to the mass matrix has 
been derived. 

Results for free vibration and buckling are in good agreement with those in 
the literature. It was also shown that for unidirectional laminated beams, warping 
cannot be ignored unless the laminate is made entirely of 0° or 90 plies. Although 
the natural frequencies may increase by almost 30% when not including warping, 
the first two vibration and buckling mode shapes remain unchanged. However, a 
more accurate function should be derived for thin-walled beam cells. 

When the structure is subject to a subtangential load, the value of the non- 
conservative parameter 77 becomes very important in design. The reason is that 
by knowing how nonconservative the load is, one can predict if the instability will 
occur due to divergence or flutter. In the case of isotropic and laminated beams 
with constant cross section, and in the absence of damping, all value of 77 < 0.5 will 
be governed by divergence instability, and by flutter instability otherwise. 


7.2 Response of Uncertain Laminated Beams 

Monte Carlo Simulation has been applied to laminated beams with randomness in 
ply orientation and the modulus of elasticity in the x-direction to study their effect 
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on the free vibration and stability of the structure. At least ten thousand realiza- 
tions of the Monte Carlo sampling have been performed to improve the accuracy of 
the analysis. A second-order Sensitivity-Based Monte Carlo Simulation (SBMCS) 
has been developed using perturbation methods. Using Taylor Series expansion, 
the eigenvalues have been expressed as probabilistic quantities. The accuracy of 
the free vibration and stability response has been compared to that obtained by 
the Exact Monte Carlo Simulation. A third approach, called the probabilistic finite 
element approach (PFEA), was also developed. It gave results that were in good 
agreement with those given by SBMCS and gave a very good prediction of the 
behavior of the fundamental natural frequency in the presence of uncertainties in 
ply angles and in Young’s modulus in the x-direction. 

The two methods employed, SBMCS and PFEA, are advantageous over simula- 
tion techniques, such as MCS, because the eigenvalue problem is solved only once. 
Also, an elegant way to obtain sensitivity derivatives was used. Based upon the 
results, the SBMCS and PFEA result in a great computational saving when one 
is interested in predicting the statistics of the fundamental natural frequency of 
laminated beams. As an example, the MCS for the case of nonconservative loading 
took about 9-10 hours whereas the other two methods took about one minute. 

For the case of free vibration, it was observed that variations in the modulus of 
elasticity in the x-direction has a greater influence over the variation in eigenfre- 
quencies when compared to the effect of uncertain ply angles. We also performed 
the reliability of beams under uncertainties in ply angles and the Young’s modulus 
in the x-direction. The reliability analysis showed that for the types of problems 
solved in this dissertation, a deterministic approach, using the traditional safety 
factor of 1.5, would have been sufficient. 
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7.3 Future Work 


Nowadays, commercial codes such as ABAQUS axe capable of analyzing the sta- 
bility of laminated structures including various types of nonlinearities. One of the 
greatest features that exists in ABAQUS is that one can employ one’s own element 
to analyze problems. It is suggested that the present formulation, and/or future 
formulations, be modified such that they can be integrated into ABAQUS. The 
drawback is that it does not have the capability to perform probabilistic analysis. 
As a future study, one could develop a computer program that would work as an 
interface to ABAQUS to analyze uncertain systems. 

The present work assumed that the type of laminated structures to be analyzed 
could be modeled using a one-dimensional model. However, real structures may 
be too complex to be modeled as beams. Thus a more rigorous study should 
expand the present study using plate and shell theory. In doing so, one can perhaps 
model and study a wide class of wings. Moreover, the use of a higher order theory 
will not require the approximations for the shear correction factors. Thus it is 
suggested that the present formulation be extended by using a higher order plate 
and/or shell theory. The present work has shown the importance of the warping 
effect in laminated structures. However, it is suggested that the structural model 
be improved by taking into account more accurate warping terms (Librescu and 
Khdeir, 1987). 

We only performed linearized stability analysis. However, post-buckling analysis 
can give us more information about the stability of the structure. The formulation 
for the post-buckling analysis has been included but not used because it was beyond 
the scope of this work. It would be most interesting to see how uncertainties affect 
the post-buckling behavior of laminated structures. 

When analyzing uncertain structures, creating a proper sampling for the purpose 
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of Monte Carlo Simulation is an extremely important step. Some researchers have 
developed new and better methods for creating these samplings than those used in 
this dissertation. Among these methods is the Latin Hypercube Sampling, widely 
used at Sandia Laboratories (Wyss and Jorgensen, 1998). Thus it is recommended 
that these new methods of generating appropriate samples be integrated into the 
analysis of uncertain systems. 

In this work, we have assumed that the laminae parameters are spatially uncor- 
related. However, the laminae parameters in real problems may be spatially cor- 
related. Spatially correlated random variables are those whose randomness varies 
from point-to-point in the structure, whereas spatially uncorrelated random vari- 
ables are those whose uncertainty is the same at any given point in the structure. 
Only a few researchers have studied the effect of having random parameters spa- 
tially correlated. Singh and Abdelnaser (1992) used a generalized modal approach 
to solve the equations of motion of a laminated composite beam, and calculated the 
random response of beams subjected to temporally and spatially correlated loads. 
Elishakoff et al. (1995) developed an exact solution to the bending of a beam with 
spatially correlated stochastic stiffness. However, no work has been found on the 
effect of spatial correlation of laminae parameters on the stability of laminated 
nonconservative systems. Thus, it is suggested that such work be pursued. 

Structures in general are made of viscoelastic materials. In fact, viscoelasticity 
can be an important factor in laminated composites (Hammerand, 1999). Moreover, 
damping may have a strong effect in the stability of laminates structures subject 
to subtangential loading. Thus it is recommended that the present analysis be 
extended to the inclusion of viscoelastic, or damping, effects. 

We assumed that the uncertain parameters can be treated as random variables 
with known (or assumed) probability distributions. However, is some cases this may 
not be true. The uncertainty, in general, is due to imprecise information and/or 
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the fact that statistical data cannot be obtained. Thus it is suggested that non- 
probabilistic approaches, such as fuzzy sets, be considered (Chen, 2000; Elishakoff 
et al., 2001). 

The work done in this dissertation can be applied to the design of laminated 
composite beams by preforming a reliability assessment of structures (Ang and 
Tang, 1984; Warren, 1997; Haidar and Mahadevan, 2000b). This will enable us to 
determine how critical the variations in the mechanical and geometrical properties 
are to the structure’s reliability. Thus it is also recommended that the design of 
laminated structures be performed using a reliability based optimization (Ba-abbad 
et ah, 2002). 



Appendix A 
Strain-Gradient Matrix 

Expressions 


In chapter 2 expressions for the Green-Lagrange strain components were given 
as 
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and were rewritten in the 

i quadratic form 
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where the displacements 

gradients in vector form are 
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The vectors h t ’s are sparse 9x1 vectors: 

h[ = {1 00000000} 

= { 000010000 } 
hi = {00000000 l} 
hj = { 000001010 } 
hi = { 001000100 } 
hg = {o 1 0 1 OOOOo} 

The matrices H,’s are very sparse 9x9 symmetric matrices: 
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Appendix B 
Laminate Constitutive Equations 


In chapter 2, the laminated constitutive equations were discussed. Here the 
coefficients for this matrix are obtained with the help of the symbolic manipula- 
tor program MATHEMATICA. Also, the coefficients used to calculate the shear 
correction factor are also presented. 

B.l Stress- Strain Relationship 

The materials considered in the present task obey Hooke’s Law: 

Sij = Cijki eia (B-l) 

where S are the Cauchy’s stress components, C^ki the 81 independent elastic 
constants, and eki the infinitesimal strain components. Symmetry due to moment 
equilibrium (S tj = S Jt ) and symmetry in the infinitesimal strains (e lJ = e ]t ) reduces 
material coefficients from 81 to 36 independent coefficients. Hyperelastic materi- 
als, for which an elastic potential or strain energy density function exists, have 
incremental work per unit volume of dW = Sjdej. This leads to 

**L.C, (B.2) 

dej dei dei dej 

Thus the material coefficients are symmetric, leading to 21 independent elastic 
constants (triclinic material). The stress-strain relationship can be further reduced 
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to 13 independent elastic constant if one plane of symmetry exists (monoclinic 
material). When two mutually orthogonal planes of material symmetry exist, the 
stress-strain relationship reduces to 9 independent elastic constants (orthotropic 
material). 

After considering a state of plane stress and condensing out e zz from the stress- 
strain relationship, Eq. (B.l) becomes 

Sij = Qijkl Zkl (B-3) 


or in indicial notation 


Si - Qij ej i,j = 1,2, 4, 5, 6 (B.4) 

Further, these stresses can be rotated about an arbitrary angle 8 and expressed in 

terms of the invariant stiffness coefficients as follows (Jones, 1999): 

Q n = Ui + U 2 cos(2 8) + U 3 cos(46>) 

Q 12 = U 4 -U 3 cos(4 8) 

_ f/ 2 sin(20) + 2 U 3 sin(40) 

<?16 = 2 

Q 22 = U 1 -U 2 cos(2 8) -I- U 3 cos (40) 

_ U 2 sin(20) - 2 U 3 sin(4^) 

Q26 = 2 

Qe 6 = - Us cos(40) 

Qu = U e + U 7 cos(28) 

Q 45 = -U 7 sin(28) 

Q 55 = U 6 -U 7 cos(28) 
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The invariant stiffness coefficients are defined as follows: 


U x = 

3 Qn + 3Q22 4- 2 Q n + 4Q66 

8 

u 2 = 

Qu - Q22 
2 

U z = 

Qn " 1 " Q22 ~~ 2 Q 12 ~ 4Qee 

8 

U 4 = 

Qll 4" Q22 + 6 Ql 2 “ 4(566 

8 

U b = 

Qll *4" Q 22 “ 2 Q 12 4* 4Q66 

8 


Q 44 + Qhh 

Ue — 

2 

U 7 = 

Q44 - Q55 


The integration of these stresses results in the extensional, extensional-bending, 
and bending stiffness matrices. These are implemented into the FORTRAN code 
as follows: 


! 


j 


! 


! 


SUBROUTINE compositel (ielm, Nlam, Exx, Eyy, nuxy, 

& nuyx, Gyz, Gxy, Gxz, ang, zk, thick, rl, r2, 

& A, B, D) 


IMPLICIT none 

DOUBLEPRECISION, PARAMETER :: zero=0.0D+0, one=1.0D+0, T0L=1.0D-7 
DOUBLEPRECISION, PARAMETER :: T0L2=1.0D-10 

INTEGER, INTENT (IN) :: ielm, Nlam 
INTEGER :: ii, jj, kk 
Material properties 

DOUBLEPRECISION, INTENT (IN) :: Exx(Nlam), Eyy (Nlam) , nuxy (Nlam) , 
t nuyx (Nlam), Gyz (Nlam), Gxy (Nlam), Gxz (Nlam) , 

t ang (Nlam) , zk(Nlam+l) , thick, rl, r2 

DOUBLEPRECISION :: C0S2ang, C0S4ang, SIN2ang, SIN4ang 
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DOUBLEPRECISION :: zkl, zk2, zk3 
! ! Laminated beam characteristics 

DOUBLEPRECISION :: Qll, Q22, Q66, Q12, Q44, Q55 
DOUBLEPRECISION :: ul, u2, u3, u4, u5, u6, u7 
DOUBLEPRECISION :: Qb(5,5), Qbbl(5,5), Qbb2(5,5) 

! ! A,B,D matrices 

DOUBLEPRECISION, INTENT (OUT) :: A(5,5), B(3,3), D(3,3) 

! 

j i ! i M M j M !! M !!!!!!!!!!!!!! ! 
j ! 

! Initializing matrices: ! 

! A, B, D, Qb ! 

! and A44.A45.A55 ! 

! ! 
j j II J J J j I j !! ! Ml M !!!!!!!!!!! ! 

Do j j— 1,5 

Do kk=j j ,5 

A ( j j ,kk) = zero 
Qb(jj.kk) = zero 
If C j j<=3 .and. kk<=3)then 
B(jj.kk) = zero 
D(jj.kk) = zero 
Endlf 
EndDo 
ElndDo 

j u !! j M ! j M II MM.!!!!!!!!!!! ! 

i ! 

! Calculation of A, B, D, ! 

! ! 

;;; i ; l l j j j! j I !!!!!!!!!!!!!!!! ! 

Do ii=l,Nlam 

i ; 1 1 ; i ; j 1 1 !!!!!! j !!!!!!!! ! 

! ! 

! Qij '• 

! ! 
n j j t MM!! mm M !!!!!!!! ! 

Qll = Exx(ii)/( one - one*nuxy(ii)*nuyx(ii) ) 

Q12 = Eyy(ii) *nuxy(ii)/( one - one*nuxy(ii) *nuyx(ii) ) 
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Q22 = Eyy(ii)/( one - one*nuxy(ii)*nuyx(ii) ) 

Q44 = Gyz(ii)*one 

Q55 = Gxz(ii)*one 

Q66 = Gxy(ii)*one 

1 1 1 11 1 11 11 t 1 11 1 i 1 1 1 1 11 11 11 


! ui ! 

! ! 

!!!!!!!!!!!!!!!!!!!!!!!!!! 

ul=( 3 . 0D0*Q11 + 3 . 0D0*Q22 + 2.0D0*Q12 + 4.0D0*Q66 )/8.0D0 
u2=( Qll - Q22 )/2.0D0 

u3=( Qll + Q22 - 2.0D0*Q12 - 4.0D0*Q66 )/8.0D0 

u4=( Qll + Q22 + 6 . 0D0*Q12 - 4.0D0*q66 )/8.0D0 

u5=( Qll + Q22 - 2.0D0*Q12 + 4.0D0*Q66 )/8.0D0 

u6=( Q44 + Q55 )/2.0D0 

u7=( Q44 - Q55 )/2.0D0 

1 1 1 1 1 11 1 l 1 1 1 1 1 l l 1 11 t 1 1 1 11 ; 


C0S2ang = dcos(2.0D0*ang(ii)) 

C0S4ang = dcos(4.0D0*ang(ii)) 

SIN2ang = dsin(2 . 0D0*ang(ii) ) 

SIMang = dsin(4.0D0*ang(ii)) 

Qb(l,l) = ul + u2*C0S2ang + u3*C0S4ang 
Qb (1,2) = u4 - u3*C0S4ang 

Qb(l , 3) = ( u2*SIN2ang + 2 . 0D0*u3*SIN4ang )/2.0DO 
Qb(2,2) = ul - u2*C0S2ang + u3*C0S4ang 
Qb(2,3) = ( u2*SIN2ang - 2 . 0D0*u3*SIN4ang )/2.0D0 
Qb(3,3) = u5 - u3*C0S4ang 
Qb(4,4) = u6 + u7*C0S2cUig 
Qb(4,5) = - u7*SIN2ang 
Qb(5,5) = u6 - u7*C0S2ang 
!!!!!!!!!!!!!!!!!!!!!!!!!! 


A, B, D 
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!!!!!!!!!!!!!!!!!!!!!!!!!! 
zkl = zk(ii+l)-zk(ii) 

zk2 = ( zk(ii+l)**2 - zk(ii)**2 )/2.0D0 
zk3 = ( zk(ii+l)**3 - zk(ii)**3 )/3.0D0 
Do jj— 1.3 

Do kk= j j , 3 

A(j j ,kk) = A(jj,kk) + qb(j j ,kk)*zkl 
B(jj,kk) = B(j j ,kk) + Qb(jj,kk)*zk2 
D(jj.kk) = D(j j ,kk) + Qb(jj ,kk)*zk3 
if C ii == Nlam)then 

if ( dabs(A(j j ,kk)) <= TOL ) A(j j ,kk)=zero 
if( dabs(B(j j ,kk)) <= TOL ) B(j j ,kk)=zero 
if ( dabs(D(j j ,kk)) <= TOL ) D(j j ,kk)=zero 
A(j j ,kk) = A(j j ,kk)/(Eyy(l)*thick) 

B(j j ,kk) = B(j j ,kk)/ (Eyy(l)*thick**2)/rl 
D(j j ,kk) = D(j j ,kk)/(Eyy(l)*thick**3)/rl**2 
A(kk, j j) = A(jj,kk) 

B(kk, j j) = B(jj.kk) 

D(kk, j j) = D(jj.kk) 

Endlf 

EndDo 

EndDo 

A(4,4) = A(4 ,4) + Qb(4,4)*zkl 
A(4,5) = A(4,5) + Qb(4,5)*zkl 
A(5 , 5) = A(5,5) + Qb(5,5)*zkl 
if Cii == Nlam)then 

if ( dabs (A (4, 4)) <= TOL ) A(4,4)=zero 
if ( dabs (A (4 ,5)) <= TOL ) A(4,5)=zero 
if ( dabs (A (5, 5)) <= TOL ) A(5,5)=zero 
A(4 , 4) =A (4 , 4) / (Eyy ( 1) *thick) 

A(4 , 5) =A (4 , 5) / (Eyy ( 1) *thick) 

A(5,5)=A(5,5)/ (Eyy(l)*thick) 

A(5,4) = A(4,5) 

Endlf 

EndDo 

return 

END SUBROUTINE composite 1 
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B.2 Equivalent Bending- Stiffness Matrix 

The bending-stiffness matrix D c is partitioned as follows: 


D c = 


Di,i Diji 
Dii.i DiihJ 

[ 0 ] 


[0] 

[D C J 


(B.5) 


where D u corresponds to the strains we are keeping, Du n the strains we are 
condensing out, and Di,n the coupling between these: 


Di,i = 


An 

A 16 

Bn 

B 16 

Ai6 

A66 

B 16 

B 66 

B n 

B\6 

Du 

Di6 

B 16 

B 66 

D 16 

D 66 


Di,ii = 


A\2 B\2 

A.26 B'ifi 

B\2 D\2 

i?26 D26 


= D 


ii, i 


Diiii = 


A22 B22 
B22 D22 


The reduced form of the first submatrix in the bending-stiffness matrix is calculated 

Dr = Dij — Diji Djjji Dii t i (B-6) 


as 
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Using the above expression, an equivalent bending-stiffness matrix D c for a 
unsymmetrically laminated beam is found: 


^Cll 

Dc\2 

Deis 

D Cl4 

0 

D C \2 


&C23 

D C 24 

0 

D C13 

&C23 

^C33 

D C34 

0 

Dcu 


-^C34 

Dc 44 

0 

0 

0 

0 

0 

K S D, 


(B.7) 


The coefficients of the above matrix are given in their analytical form. Following 
is the subroutine included in the FORTRAN program NLbeam21.f b 


SUBROUTINE DcMatrix(A, B, D, rl, r2, ShearFac, Dc) 
IMPLICIT none 

DOUBLEPRECISION, PARAMETER :: zero=0.0D+0, one=1.0D+0, 

& T0L=1 . OD-7 

INTEGER :: ii, jj, kk 

DOUBLEPRECISION, INTENT (IN) :: ShearFac, A(5,5), B(3,3), 
& D(3,3) , rl, r2 

DOUBLEPRECISION, INTENT (OUT) :: Dc(5,5) 

DOUBLEPRECISION :: delta 


Initializing matrix Dc 


delta = B(2,2)**2 - A(2,2)*D(2,2) 

Dc (1 , 1) = ( A(2,2)*B(1,2)**2 - 
& 2 . 0D0*A(1 ,2)*B(1 , 2) *B(2,2) + 

& A(1,1)*B(2,2)**2 + A(1 ,2)**2*D(2,2) - 


1 Note that all “3” indices in the A, B, D matrices represent “6”, i.e., -4(1 , 3) -> 4 16 . 
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ft 

ft 

E 

ft 


Dc(l ,4) 


Dc(l , 5) 
Dc(2,2) 


t A(l, 1)*A(2,2)*D(2,2) 

U ) /delta 

Dc(l ,2) = ( -A(2,3)*B(1,2)*B(2,2) + A(1,3)*B(2,2)**2 + 
fc A(2,2)*B(1,2)*B(2,3) - A(1,2)*B(2,2)*B(2,3) 

fe A(1,3)*A(2,2)*D(2,2) + A(1,2)*A(2,3)*D(2,2) 

ft ) /delta 

Dc(l ,3) = ( -B(1,2)**2*B(2,2) + B(1,1)*B(2,2)**2 + 
ft A(2,2)*B(1,2)*D(1,2) - A(1,2)*B(2,2)*D(1,2) 

A(2,2)*B(1,1)*D(2,2) + A(1,2)*B(1,2)*D(2,2) 

) /delta 

= ( B(1,3)*B(2,2)**2 - B(1,2)*B(2,2)*B(2,3) - 
A(2,2)*B(1,3)*D(2,2) + A(1,2)*B(2,3)*D(2,2) 

A(2,2)*B(1,2)*D(2,3) - A(1,2)*B(2,2)*D(2,3) 

) /delta 
= zero 

= ( A(3,3)*B(2,2)**2 - 
2 . 0D0*A (2 , 3) *B (2 , 2) *B(2 , 3) + 
t A(2,2)*B(2,3)**2 + A(2,3)**2*D(2,2) - 

t A(2,2)*A(3,3)*D(2,2) 

t ) /delta 

Dc (2 , 3) = ( B(1,3)*B(2,2)**2 - B(1,2)*B(2,2)*B(2,3) 
k A(2,3)*B(2,2)*D(1 ,2) + A(2,2)*B(2,3)*D(1 ,2) 

ft A(2,3)*B(1 ,2)*D(2,2) - A(2,2)*B(1 ,3)*D(2,2) 

ft ) /delta 

= ( -B(2,2) *B(2,3) **2 + B(2,2)**2*B(3,3) + 
A(2,3)*B(2,3)*D(2,2) - A(2,2)*B(3,3)*D(2,2) 

A(2,3)*B(2,2)*D(2,3) + A(2,2)*B(2,3)*D(2,3) 

) /delta 
= zero 

= ( B(2,2)**2*D(1,1) - 
2 . 0D0*B ( 1 , 2) *B (2 , 2) *D ( 1 , 2) + 

& A(2,2)*D(1,2)**2 + B(1 ,2) **2*D(2,2) - 

& A(2,2)*D(1.1)*D(2,2) 

ft ) /delta 

Dc(3,4) = ( -B(2,2)*B(2,3)*D(1 ,2) + B(2,2)**2*D(1,3) 

& B(1,2)*B(2,3)*D(2,2) - A(2,2)*D(1 ,3)*D(2,2) 

& B(1 , 2)*B(2 ,2) *D(2,3) + A(2,2)*D(1,2)*D(2,3) 

& ) /delta 


Dc(2,4) 


ft 


Dc(2, 5) 
Dc(3,3) 
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Dc(3,5) = zero 

Dc(4,4) = ( B(2,3)**2*D(2,2) - 
k 2 . 0D0*B(2,2)*B(2,3)*D(2,3) + 

k A(2,2)*D(2,3)**2 + B(2,2)**2*D(3,3) - 

k A(2,2)*D(2,2)*D(3,3) 

k )/delta 

Dc(4,5) = zero 

Dc(5,5) = ShearFac*( A(5,5)- A(4,5)**2/A(4,4)) 

Do ii=l,5 

Do j j=ii,5 

if ( dabs(Dc(ii , j j ) ) <= TOL ) Dc(ii, j j)=zero 
Dc(jj ,ii)=Dc(ii, j j) 

EndDo 

EndDo 

! 

return 

END SUBROUTINE DcMatrix 


B.3 Transverse Shear Coefficient Factor 


In order to calculate the transverse shear strain energy density, we need to calculate 
the following matrices and vectors laminawise. Thus here we present the final results 
obtained using MATHEMATICA. 


The A is a 4 x 4 matrix defined as 


A = 


An 

-A- 12 


A12 


A22 


Au(l,l) 

0 

Ai 2 (1,1) 

A ia (l,2) 

0 

0 

0 

0 

Aia(l, 1) 

0 

A-22(l) 1) 

A.22(E 2) 

A 12 (l, 2) 

0 

A.22(l) 2) 

A 22 (2,2) 


(B.8) 


Let N\ &m be the total number of plies considered. Then the 2x2 matrix Ay can 
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be defined as 


V2 

f <*j 

Ly= J ~G- 

—h/2 


dz 


(B.9) 




■^lam 

A 12 (1,1) = E 

k = 1 



4 + i-4\ B k n D k 33 


;\ g i 
J 41 


+ 


+ 


+ 


4 4 > 

2 fc+l - Z k 


z k + 1 “ 


2 2 ' 
2fc+l - Z k 


4G* 

D& D& f\k 

&11 n 33 ^11 ^33 


2G k 


2G k 


Mi Bj 3 _ ^4*3 _ ^ 4a h 2 B k n P k 33 

yk 4G fc 4G k 8G k 


G* 


hA k n B k 3 , h 2 4i 4 fc 3 , h 2 A k n D k . 

I _ “T 


33 


G* 


8 G 


8 G 


+ (2/c+i — 2:*;) ( — 


h 2 A k n B k 33 fi 3 ^B 3 fc 3 , ^4*3 _ 

4G* 16 G* 16G fc 64 G k 
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^lam 

Aw( 1,2) = £ 

/c=l 


+ 


+ 


+ 


-k + 1 


- Zl. 


z k+i 


- Zl. 


z k+l 


- Zl 


r 3 - 7 3 ' 

z k + 1 ffc 


+ (2fc+l — Zfc) 


B k n Dv 


13 



B k n B \ 3 , ^>13 



' h 2 _ ft 3 Bfl ^13 _ ^llffil + ^ fll ft 


4G fc 


16 G 


16 G 


64 G 



•^lam 

■A.22(l) 1) = 

fe=l 

+ (^ 


4«-4\ 



33 


4G 


G 



h 2 D k 33 D k M 


8G 


h 4 -D33 G3 3 

64G fc 
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■^lam 

a 2 2(i, 2 ) = y: 

fc=l 


+ 


+ 


+ 



4 4 s 

z k + 1 - 


4+i - 4' 


2 2 ' 
2fc+l - z k 


4+i -4\ 

^13 ^33 

5 J 

4G fc 

*31 

<Q 

<CQ 

pfc r>fc 
^13 ^33 


2 G k 


2G k 


+ (Zfc+l — z k) ~ 


B k 13 B k 3 3 /lB 3 fe 3^13 ^3^33 

G k 4G k 4 G k 8G k / 

h B k 3 B 33 h 2 B 33 P k 3 h 2 B k 3 -P33 \ 

G fc 8G fc 8G fc y 

h 2 B k 3 B%z h^pl 3 

,fc + 16G fc 16 G* 64 G fc 


4G 


•^lam 

A 22 (2,2) = £ 

/c=l 


+ (^fc+l “ Z k) 


4+1 - 4 




B k 3 B k 3 h B k 3 P k 3 _ h 2 P k l3 P l3 
2 G k 8 G k 

hB k 3 B k 3 h 2 B k 3 P k 3 \ 

4G k 

h 2 B k 3 B k 3 fi 3 B& Dj 3 + 


8G 


64 G 


The B is a 4 x 1 vector defined as 


{£} 


= < 


' Bi(l) ' 
0 

B 2 (l) 


(B.10) 
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Then the 2x1 vector Bj can be defined as 

h/2 


Bi 


-/ 

-h/2 


a? (3 

G 


dz 


(B.ll) 



+ 


+ 


+ 


gn Dn 
4G fc 





8G 


h 4 B k n D k v 
6AG k 



B " k fyk 
11 ^33 



D k n D k 


33 


4G 

gu 

2G fc y 

hB k 3 D k u 

4G k 


h -Bji 

4G fc 


+ 


ft 2 gn g| 

8G k 


+ (^fc+l ” 2 fc) i 


h B‘, BSs . ) 

k 8 G k 8G k ) 

h 2 B k n B k 3 /r 3 4 fc 3 ^ii , ^ gn ^3 _ 

+ 16 G k 16 G k 64 G k 
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f)k r\k 
U U U \1 


+ 


+ 


+ 


+ ( z k + 1 — z k) 


AG K 



B k n Dy 


13 


2 G k ) 

hB k 3 D k u , hB k n D 

AG k 

h 2 B k 3 D k n 


u D k n D 


'13 


AG 


8 G k 


h 2 #3 


8 G k 8 G k 

h 2 B k n B k 3 h 3 B k a D k n h 3 B k n D k : 


13 




4 G k 


16 G k 


16 G* 


64 G k 


The C is a 1 x 1 matrix defined as 


C 


h / 2 

-I 

—h/2 


(3 t (3 


dz 


(B.12) 



+ 


+ 


0fi 0fi 


AG k 


+ 


B* n D* u 


V 


+ { z k + 1 — z k ) 



h* 0n 0ii 

8G k 


h 2 B k n B h n 


AG k 


2G k 

h 2 B k n D k u 
AG k 

8G fc 64 G fc 


Appendix C 
Tangent Stiffness Matrix 


In chapter 3, the tangent stiffness matrix was derived and expressed using ma- 
trices and vectors. However, by using the indicial notation, the stiffness matrices 
are easier to implement in computer codes. 

In indicial notation, the tensor or matrix components are explicitly specified. 
Thus a vector, which is a first-order tensor, is denoted in indicial notation by x l , 
where the range of the index is the number of dimensions of x , n S( i. Indices repeated 
twice in a term are summed, in conformance with the rules of Einstein notation. In 
fact, the indicial notation is almost unavoidable in the implementation of FEM; for 
programming the finite element equations, the indices must be specified (Belytschko 

et al., 2000). 

The virtual work done by internal forces can be written as 

SW M = JJJs,Se,d f (C.l) 

f 

where Sj are the dimensionless PK2 stresses and are energetically conjugate to e 3 , 
which are the dimensionless Green-Lagrange strains. 

Since the virtual displacements are defined in the same space of functions as 
the finite element space of functions, the total virtual generalized strains can be 
expressed in terms of the virtual dimensionless linear and nonlinear strains, in 
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indicial notation, as follows: 

S ei = p-Sq. = B;fH (C.2) 

dqi 

= (Bj + B% L ) S Qi (C.3) 

j = 1,2, , ^strains f = 1, 2, . . . , Nfof 

where Strains is the total number of generalized strains (iV strai ns = 5), iVdof is the 
total number of degrees of freedom of the element (JV dof = 21), B*f the strain- 
displacement matrix, B£ the dimensionless strain-displacement matrix independent 
on the displacements, and B* L the dimensionless nonlinear strain-displacement, 
matrix linearly dependent on the displacements. 

Expressing Eq. (C.l) in terms of the generalized strain and stress vectors, we 
get 

= JJ Tj 5£j dCl j = 1,2,..., JV strains (C.4) 

n 

where N Btrains is the length of the generalized strain vector and in this research it 
is equal to five (iV 8trains = 5), 8sj are the dimensionless virtual generalized strains 
given by Eq. (C.2), and T 0 the generalized stresses given by Eq. (2.55). Using Eq. 
(C.2), the internal virtual work becomes 

<5W int = Sqijj ^ T J dCl (C 5) 

^ ‘ ' 

ir 


Now, the tangent stiffness matrix is given by 

dq k JJ dq t dq k J J dq> dqk_ 


Ktk = 
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where k = 1,2,. 
integral, 


. , N do f . Using the chain rule for the underlined term in the second 


dr. 

= (&) 

de r 

= D c = D c .B? k 

(C.6) 

dq k 

\de r ) 

dq k 

dq k 


T 

1 J 

— D Cjr E r 



(C.7) 


Thus, the tangent stiffness matrix becomes 

* = irs^in 


tD C]r B? k L dCl 


I! 


+ II D Cjr B L rk dV + JJ B% L D Cjr B? k L dCl 


(C.8) 


where 


dB*f dBfi dB" L dB% L 




J 1 


(C.9) 


dqk dq k dq k dq k 
= 0 

and is independent of displacements! 

Thus, decomposition of the element tangent stiffness matrix can be written as 


Kt k = K% + Kg + K\ 


ik 


(C.10) 


where , Kg, and Kg denote the element linear, initial-displacement, and geo- 
metric stiffness matrices, respectively, where the above matrices are given by 


fl rl/2 

K% = / B^ D c . r Bg k dy dx (C.ll) 

J0 J- 1/2 

y*l rl/2 

Kg = if {2 Bg D Cjr Bg k L + B" l D Cjr Bg k L } dy dx 
Jo J- 1/2 


(C. 12) 
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where D Cjr is the dimensionless equivalent bending-stiffness matrix. 

Moreover, when non conservative forces are present we have to add the loading 
stiffness matrix to the global tangent stiffness matrix. This loading matrix can be 
defined as 


Ktk = 


djr 

9qt k 


(C.14a) 


For the case of a beam subject to a tangential follower load, 


K, 


ik 


= Nlh^ + Nl, 


dq, 


't k 


902 , , 
9qt k 


Qt r + ^Lji 02 


'jfc 


(C.15a) 


Appendix D 
Dynamic Condensation Technique 


Reduction methods have been studied and developed in the past decades. Noor 
(1994) gave an extensive review on recent advances and application of reduction 
methods. 

Here we focus on the dynamic condensation technique, a technique not com- 
monly used among researchers. The reason is that most of the problems of interest 
are static and not dynamic. However, here we highlight one method of deriving the 
dynamic condensation technique. 

In order to proceed, the displacement vector is rearranged by partitioning the 
relevant terms corresponding to external and internal degrees of freedom as follows: 

q tl = {qi,q2} r ( D1 ) 

where qi are the exterior nodes to be kept and q 2 are the interior nodes to be 
condensed. Now the load vector can be partitioned in a similar way: 

F = {Fi,F 2 } t (D.2) 

The element stiffness matrix and mass matrix are also rearranged and partitioned 
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according to Eq. (3.42) as follows: 



K„ 

k 12 

K 2 i 

K 22 

M n 

M 12 

M 21 

M^22 


The following definitions are introduced: 

qi = q 

q2 = R q + € 


(D.3) 

(D.4) 


(D.5a) 

(D.5b) 


where q is the displacement vector with only those displacements not to be con- 
densed, e is the perturbation, and R is the matrix which allows the condensation 
to take place. In order to find matrix R, Eq. (D.5) is substituted into the strain 
energy 


2Wi„ t = q T [Kn] q 4- € T Ka 2 € + (D.6) 

q T [R t K 2 i 4- K 12 R + R t K22R] q 4- 
26 t [K21 4- K22R] q 

Now, the underlined term in Eq. (D.6) must vanish to decouple the variables q and 
e for static problems. This is the basis for static condensation. This leads to an 
expression for the matrix R: 


R = — K 22 1 K2i 


(D.7) 


Since e is very small perturbation, its quadratic term in Eq. (D.6) is neglected. 
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Further, substituting Eq. (D.7) into Eq. (D.6), we get 

2Wj n i = q r [Kn 4 - R t K 2 i + K12R + R r K22R-] Q 
= q T [Kn - K 12 K22" 1 K 2 i] q 

+q T [— K 2 i T K22 - ' r K2i + K2i r K 2 2 -T K22K22 X K2i] q 
= q T [Kn - K21 T K 2 2“ T K21 - K 12 K22 _1 K2, + K 12 K22" 1 K2l] q 
= q T [Kn - K 12 K 2 2“ 1 K2i] q 

Thus, 

K R = K n - K 12 K22 _1 K2i (D.8) 


The dynamic condensation is performed by first defining the following matrix: 


Q = 


— K 22 X K21 


(D.9) 


where I is the identity matrix of same dimensions as R. 

The reduced matrices are found then by using the following relationship: 


k r 

= q t k q 

(D.10) 

f r 

= q t f 

(Dll) 

m r 

= q t mq 

(D. 12 ) 
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For example, the stiffness matrix is 

I 

-Kia^Kai _ 

Kn — K12K22 x K 2 i 
K21 — K22K22 1 K'2l 
= K11 — K12K22 ! K21 — K12K22 L I^21 + K12K22 K21 

= Kn — K12K22 ! K 2 i 

which gives the same as Eq. (D.8). A similar procedure can be followed for the 
mass matrix and load vector matrix. 


— K12K22 


-1 


k r = 


-K 2 i T K 2 2- T 


K n K 12 

K21 K22 


Appendix E 
Probabilistic Formulation 


The approach involving the perturbation method was recently used by Zhang 
and Ellingwood (1995). Haidar and Mahadevan (2000b) provide a good understand- 
ing regarding stochastic finite element analysis (SFEA). Basically, the stochastic 
finite element analysis relies on the fact that the sensitivity derivatives are known. 

The present formulation for the probabilistic finite element analysis is similar to 
that of SFEA. This usually leads to recursive equations that are solved sequentially 
to obtain the variation of eigenfrequencies, eigenvectors, and buckling loads. 

Eigenvalue derivatives 

Let K be the probabilistic linear stiffness matrix, M the probabilistic mass matrix, 
L the probabilistic loading matrix, A*, the probabilistic eigenfrequency of the k 
mode, P k the probabilistic buckling load of the k th mode, ij> k the probabilistic left 
eigenvector of the fc th mode, and 4> k the probabilistic right eigenvector of the k tb 
mode. Thus the probabilistic eigenvalue problem is expressed as 

W k } T [K-P k L-\ k M} = 0 (E.l) 

[K-PkL-Xk M] {<f> k } = 0 (E-2) 


As mentioned in chapter 5, all quantities in Eq. (E.2) are perturbed about their 
mean value. Thus, they are expressed in terms of their mean-centered zeroth-, 
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first-, and second-order rates of change with respect to the random variables as 
described in section 5 . 2 . 2 . Further let the first order rate of change of the k th mode 
right eigenvector with respect to the random variable evaluated about the mean 

be 

d{<f> fc} 


dxi 


_ - W) 


(E. 3 ) 


x=x 


and the second order rate of change of the fc th mode right eigenvector with respect 
to the z th and j th random variables evaluated about the mean be 


a 2 M 

dxi dxj 



(E. 4 ) 


Using the same notation for the eigenfrequencies and buckling loads, we can express 
the following: 


K(r lt r 2 , • ■ 

• , r n ) = K° + ^2 K i €i + 9 ti€j 

»= 1 »=i j = 1 

(E. 5 ) 

M(ri,r 2 ,.. 

. , r n ) = M° + Y, M \ * + lib it, ^ 

«=1 i=l j = 1 

(E.6) 

L(r i,r 2) . • 

■ , r n ) = L° + ^ L{ €i + eitj 

«= 1 «=i j=i 

(E. 7 ) 

A fc (ri,r 2 ,. 

• , r n ) = A° + ^ A^ Afc-J 

»=i «=i j=i 

(E.8) 

{Mn,r 2,-- 

• ,r„)} = {^fe} + 2^-' 

i= 1 *=1 J=1 

(E. 9 ) 

Pk{ri,r 2 ,• 

. r •„) = Pj + £ PL e, + 5 £ E e “i 

1 t=l j = l 

(E. 10 ) 


After substituting these perturbations into Eq. (E. 2 ), a probabilistic eigenvalue 
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problem is formulated (for simplicity, the perturbed eigenvectors were not substi- 
tuted): 

k°+j 2 k > ^ + Rl i ^ I W = 


i=l 


i— 1 j=l 


V i=i i=i j=i / 

M °+E M ! f i+^EE M " e ‘ e J W 

t=l t=l j — 1 J 

' t=l i=l j = l / 


(E.ll) 


t=l 


1=1 j = 1 


Further, we expand the above and equate all terms of f , in the expansion. The 
uncertainties are in general small and, as a consequence, in the applied perturbation 
technique it is sufficient to only consider up to second order terms. This leads to 


(E. 12) 


+ 


+ 


M 


[K° - P° k L° - A °M°] {0 fc } 

J2 [K\ - P° k Lj - PLL 0 - A LM° - A ®Mf] f-i 
. 1=1 

i=l j=l -1 

lY,i2 [~ x lii M ° - x ii M l - A t M « - A ° M */] W = 0 


i=l j=l 


Now substituting the perturbed eigenvectors, and keeping only second-order 
terms, we get 


[K° - P° k L° - \° k M°] {*»} 


APPENDIX E. 


226 


+ £ [K° - P° k L° - X° k M°] { 0 L> * 

i— 1 

+ £ [K\ - PiU - PLL° - A',M» - A°M'] Me, 
1=1 

+ E E [*' - - pI ^ L ° ~ A * M ° _ A * M <] 

t=i j=i 


+ P °* L ° - A ° M °] e ‘ e > 

1=1 j=l 

- P^Lj - P{jL{ - P° k L{; - \' ki Mj - AljMj - A °M?/] {<$} ^ 


We want the above expression to be valid for all small e. Thus we set the 
coefficients of e N (N = 0, 1,2) to zero. Moreover, for the type of problem solved 
throughout this dissertation, the loading and mass matrices do not depend on the 
mechanical properties of the beam; as a consequence, their sensitivity derivatives 
vanish. This leads to three recursive equations: 

e° : [K° - P° k L° - A°M°] {</>£} = 0 (E. 13) 

e 1 : [K° - P° k L° - A^M 0 ] {<^} (E.14) 

+ [Kl - pLl° - aLm°] {*£} = 0 

: [K\ - P' t ,L° - A',M»] {<*>',} (E.15) 

+ 1 [K° - P° k L° - A JM°] {<} 

II oil rO \II Ji/fOl r Afi 


APPENDIX E. 


22 7 


Now premultiply Eqs. (E.14) and (E.15) by {ip 2} T to get 

t 1 ; MV [K° - PlL° - AjM 0 ] {^J (E.16) 

+ MV [K‘ - pL,l° - ALM°] {K} = 0 

c 2 : 2 {i l>l) T [K\ - Pl,L" - A^M 0 ] {<) (E.17) 

+ MV [X° - P 0 t L a - AjM 0 ] «} 

+ MV [K'j - PlijL 0 - A" M“] {^} = 0 

Note that the present deterministic eigenvalue problem can be separated into two 
eigenvalue problems: 

MY [K° - P° k L° - A°M°] = 0 
[K° - P° k L° - A°M°] {<f>° k } = 0 

where the loading matrix may be unsymmetrical. Since neither the left nor the 


right eigenvectors 
equations become 

are zero, [ K° - P° k L° - A ° k M°] = 0. Thus the three recursive 

e° 

: [K° - PlL 0 - A°M°] {$} = 0 

(E.18) 

e 1 

: MY [*■; - pLl° - M) = 0 

(E.19) 

e 2 

: MY [ K{ ! - Pl!jL° - KijM 0 } MY = 

(E.20) 


- {^} T [*; - pLl" - ALM°] my 

- MY [j i' - PljL 0 - M) 
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Eigenvector derivatives 


The first-order variations in the eigenvectors can be obtained by using the method 
described by Fox and Kapoor (1968). Murthy (1986) developed several methods 
based on the generalized Rayleigh quotient for the sensitivity analysis of the eigen- 
value problem. Later, Bergen and Kapania (1988), and Kapania et al. (1991), 
presented a method for calculating the shape sensitivity of a wing aeroelastic re- 
sponse with respect to changes in geometric shape. However, these methods are 
restricted to buckling and/or eigenfrequencies. Plaut and Huseyin (1973), and re- 
cently Adhikari and Friswell (2001), have studied the derivatives of the eigenvalue 
problems for nonconservative systems. However, here we extend the method de- 
scribed by Fox and Kapoor (1968) to nonconservative systems and to the stability 
analysis using the dynamic criterion. 


We calculate the left and right eigenvector sensitivities separately and use the 
fact that the sensitivity derivatives of the mass and loading matrices are zero. 
Since the right eigenvectors form a complete set of vectors, an eigenvector can be 
represented by the linear combination of all other right eigenvectors. Thus the 
derivative of the k th mode eigenvector with respect to the i th random variable 
evaluated about the mean is represented as follows: 


d{4> fc } 

dxi 


_ o = {flw} M 


j-l 


(E.21) 


where {<pj} is the eigenvector corresponding to the j th mode, and n corresponds to 
the total number of modes (dimensions of the stiffness matrix). Thus the problem 
reduces to calculating the coefficients ajy . 


Let us start with the following eigenvalue problem: 


[K — Pfc L — Afc M] {tf> k } = 0 


(E.22) 
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Differentiating the above equation with respect to Xi, a random variable, using 
Eq. (E.21), and rearranging the terms, we get 

[K - P k L - A* M] £ ag {4> j } = ~ [K{ - P f ki L - M] {<M (E.23) 

j = l 

Premultiplying the above equation by the transpose of the left eigenvector, {^ m } , 
we get 

E “8 (W T |K-P*1- a* M] {</.,-} 

j=i (E.24) 

= - {0„} T [Kl -PL L - At M] {<M 

To normalize the eigenvectors we need two independent criteria. Thus let us nor- 
malize the eigenvectors such that 

{V’j} [A'l] {<^j} 1 and { } n th nonzero element nonzero element 

for a selected value of n. As a result, 

{0,} T (K-ftL]{0 j } = A, 

Moreover, for distinct eigenfrequencies, the right and left eigenvectors satisfy biorthog- 
onality criteria, i.e., 

{0„,} T |M] {<*>,} = 0 {V'm} 7 [K - ft L) = 0 Vm#j 

Thus Eq. (E.24) becomes 

E4> {$j} T [K — Pfc L — AfcMj {0,} 

j - 1 

= -{yJ T [a:'-pLl-aLm]{0 1 } 


(E.25) 
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Using the biorthogonality criteria, the constants ajy ’ s for all j k are found as 

_(0 _ {*,} T [Kl - PW {M (E.26) 

“« " A* - A, 

To find the constant af$% we use the two normalization criteria defined above. Let 
us take the derivative of the first normalization criteria with respect to a random 
variable x t \ 

[M] {4>,} + {* l f [M'] {<*>,}_+ {+,}* [M] {4>L} = o (E.27) 

s V 

=0 

Since the left eigenvectors form a complete set of vectors, the eigenvector can be 
represented by a linear combination of all other eigenvectors: 

{^L} = S 6 ${^} (E - 28) 

j = i 

Substituting Eqs. (E.21) and (E.28) into Eq. (E.26), we get 

a jj + b n = 0 (E.29) 

If the first element of the left and right eigenvector remain equal, then so do the 
corresponding elements of the derivatives. Thus by differentiating the second nor- 
malization criterion and using Eqs. (E.21) and (E.28), we get (Adhikari and Friswell, 
2001): 

aji - bjj = 0 (E.30) 

Thus, ajj = bjj = 0. 
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